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I. Introduction 

This final report will attempt to concisely summarize the activities and 
accomplishments associated with NASA Grant NAG-1 -793. The project started on 
July 1, 1987 and officially terminated on December 31, 1994. While the total 
funding for the project was $1 10,395, many lengthy periods existed in which little or 
no funds were available for expenditure by the project; and all grant funds were 
essentially expended by August 31, 1993. Fortunately, the effort was maintained by 
significant financial support by the Aerospace Engineering Department in the form 
of Graduate Assistantship funds and faculty salary support and by moral and 
technical support from NASA Langley. In spite of these difficulties, significant 
accomplishments were achieved by the project; and these are summarized below. 

II. Personnel 

The individuals who have been associated with the project are as follows: 

Leland A. Carlson, Professor of Aerospace Engineering — Dr. Carlson 
served as the principal investigator for the project. At various times, Dr., Carlson 
was partially supported by the project. 

Hesham M. El-banna, Graduate Research Assistant and Graduate Assistant 
Non-Teaching (GANT) - Hesham El-banna joined the project at its inception. 
During the project, he earned his Masters' and Ph.D. degrees using research 
associated with the project for his thesis and dissertation. Dr. El-banna was 
partially supported by the project at various times. He was also extensively 
supported by the Aerospace Engineering Department as a GANT. 

Alan Arslan - Graduate Research Assistant Non-Teaching -- Alan Arslan 
joined the project in Fall 1992. He was supported by the Aerospace Engineering 
Department as a GANT. He used research associated with the project for his 
masters' thesis. 


III. Accomplishments 

The primary accomplishments of the project are as follows: 

1. Using the transonic small perturbation equation as a flowfield model, the 
project demonstrated that the quasi-analytical method could be used to obtain 
aerodynamic sensitivity coefficients for airfoils at subsonic, transonic, and 
supersonic conditions for design variables such as Mach number, airfoil thickness, 
maximum camber, angle of attack, and location of maximum camber. The approach 
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was validated by comparison to results obtained using the finite difference 
technique. It was established that the quasi-analytical approach was an accurate 
method for obtaining aerodynamic sensitivity derivatives for airfoils at transonic 
conditions and usually more efficient than the finite difference approach. These 
initial results were among the first aerodynamic sensitivity results obtained by the 
quasi-analytical approach for transonic conditions. 

2. The usage of symbolic manipulation software to determine the appropriate 
expressions and computer coding associated with the quasi-analytical method for 
sensitivity derivatives was investigated. Using the three dimensional fully 
conservative full potential flowfield model, it was determined that symbolic 
manipulation along with a chain rule approach was extremely useful in developing a 
combined flowfield and quasi-analytical sensitivity derivative' code capable of 
considering a large number of realistic design va riables . Various methods of 
solving the resulting large system of quasi-analytical equations were investigated. It 
was concluded that for the direct solver approach, that the iterative conjugate 
gradient method was accurate, capable of handling a large number of design 
variables, and more efficient than the finite difference approach. 

3. Using the three dimensional fully conservative full potential flowfield model, 
the quasi-analytical method was applied to swept wings (i.e. three dimensional) at 
transonic flow conditions. The study included as basic design variables freestream 
Mach number, wing angle of attack, airfoil thickness, airfoil camber, location of 
airfoil maximum camber, wing twist angles at four spanwise locations, and the 
coordinates of the wing tip. From sensitivity derivatives for these design variables, 
sensitivities were also obtained for other variables of interest such as wing area, 
aspect ratio, wing sweep, and taper ratio. The resultant sensitivity derivative results 
were verified by comparison with finite difference computations. Sensitivity 
derivatives were obtained chordwise for 3 C p /3Xd, spanwise for 9Cj/3Xp , and overall 
for dC l /3X d , where Xq is any design variable. The sensitivity derivatives were also 
use to predict pressure distributions and aerodynamic coefficients at conditions 
different from those at which the derivatives were obtained. These predicted results 
demonstrated that sensitivity derivatives could be used over limited ranges for 
predictive purposes. Sensitivity derivatives were also obtained over a range of 
Mach numbers ranging from 0.8 up to 1.2 and for a variety of wing airfoil sections. 
These results demonstrated the feasibility and usefulness of the quasi-analytical 
approach for obtaining aerodynamic sensitivity derivatives about wings. It is 
believed that they were among the first sensitivity results to be obtained using the 
quasi-analytical method for wings at transonic conditions. 

4. The incremental iterative technique has been applied to the three 
dimensional transonic nonlinear small perturbation flowfield formulation, an 
equivalent plate deflection model, and the associated aerodynamic and structural 
discipline sensitivity equations; and coupled aeroelastic results for an aspect ratio 
three wing in transonic flow have been obtained. This approach permitted the use 
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of finer grids and inclusion of both aerodynamic and structural sensitivity derivatives 
with the final results including full aeroelastic coupling. In addition, system 
sensitivity derivatives were obtained. Results were obtained for nine aerodynamic 
design variables and four structural design variables. The results demonstrated the 
usefulness and feasibility of combining the incremental iterative approach with the 
quasi-analytical formulation for obtaining both discipline and system sensitivity 
derivatives. Again, these are among the first results to utilize the quasi-analytical 
approach to obtain aerodynamic-structural coupled sensitivity derivatives and 
system sensitivity derivatives. However, it appears that further studies are needed 
in methods associated with determining system sensitivities and of utilizing this 
information in optimization procedures. This effort is discussed in Section IV below. 

IV. Progress in the Last Six Months 


During the past six months, Arslan and Carlson as part of a pilot study have 
applied the incremental iterative technique to the transonic nonlinear small 
perturbation formulation, an equivalent plate deflection model, and the associated 
discipline sensitivity equations, to obtain coupled aeroelastic results for an AR=3 
wing in transonic flow. This integrated approach allows the use of finer grids and 
simultaneously yields the aerodynamic and structural deflection solutions, the 
aerodynamic sensitivity derivatives for nine aerodynamic design variables, the 
structural sensitivities for four design variables, and the coupling derivatives needed 
for the system derivatives, which are computed subsequently. It is outlined in Fig. 1 . 


[Sweep flowfield by cross planes] 

[Solve flowfield equations for A<t>jj k ] 

j^Solve aerodynamic sensitivity eqs. for A^^j 

[Solve structural eqs.for A5ij k ] 

Solve structural sensitivity eqs. for A^g|^J 

Solve coupling derivative eqs. for A ^^ - - j j 

Solve coupling derivative eqs. for A ^ 3AC a5 ^ J J 

[Update aerodynamic and structural boundary conditions] 
[Iterate until convergence] 


Solve system sensitivity eqs. for etc. 


Fig. 1 -- Integrated Solution Approach 

Basic aerodynamic sensitivities were obtained at all 97x16x16 flowfield points 
while aerodynamic coefficient and structural deflection derivatives were computed 
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at 980 wing surface points, comprising ten half-span stations each having 49 points 
on top and 49 on the lower surface. System derivatives were obtained with the 
Global System Equation Method 2 (GSE2). Since Sobieski has suggested that 
system derivatives can be obtained using condensed information, system 
sensitivities were computed from the fine grid aerodynamic, structures, and 
discipline sensitivity resul ts for two condense d cases. The first computed system 
sensitivities at eight half-span stations each with 13 chordwise locations and the 
second used five half-span stations each with 25 chordwise values. Limiting the 
number of system sensitivities is desirable since the coupling derivatives required to 
compute them treat each 5 and AC p considered in the system formulation as a 
design variable. Thus, even for the condensed problem using a 25x5 system grid, 
the number of design variables was effectively 138; and 125 lengthy coupling 
vectors a5/aAC pk and 8AC p /a5 k had to be computed. Obviously, the inclusion of 
system sensitivities gre atly i ncreases the problem complexity. While results were 
obtained at each wing statio n, representative results for the 97x16x16 flowfield, 
49x10 structural, and 25x5 system case are shown on Fig. 2. The wing was at 

= 0.82, a = 2°, had 1° of twist (T, ip ), an d the airfo il s ections varied f rom a NACA 
2406 at the root to a NACA 1706 at the tip. The deflected wing position is shown 
dotted, the dC p /3Tti P curves are for the upper and lower surfaces, and in the two 
lower plots the discipline derivatives are dotted while the system sensitivities are 
solid. Note that the flow is supercritical, that the wing has significant deflection and 
twist due to aerodynamic loading, and the upper surface shock wave strongly 
affects the aerodynamic derivatives. Also, note that the system derivative dAC p /dT tip 
is significantly different from the discipline value near the trailing edge due to twist 
induced by aerodynamic-structural coupling, and that d5/dt, where t is a wing 
structural thickness parameter, differs in magnitude and sign from the 
corresponding discipline result, 85/5t. 

Unfortunately, comparison of the 13x8 and 25x5 system derivative results 
indicates differences in values, magnitudes, and sometimes signs. Since these 
sensitivities were obtained from the same fine grid aerodynamic, structural, and 
discipline sensitivity solutions and since the differences do not appear to be due to 
numerical error, they must be associated with the system derivative solution 
approach, the number and location of condensed points considered, etc. 

Since accurate system derivatives are required before the optimization portion of 
a MDDO process can be applied to transonic wing design, the methodology and 
approach for computing system derivatives for a transonic aeroelastic wing needs to 
be further investigated. In addition, since the present study utilized a pilot code and 
was primary a research investigation of feasibility, further work is needed to develop 
an aerodynamics flowfield solver and sensitivity module that is suitable for 
engineering applications and studies. 

A copy of Mr. Arslan's masters' thesis and an abstract of a proposed AIAA 
paper, which discuss this effort and include many of the details, will be sent under 
separate cover to the project monitor. 
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V. Publications and Degrees 

The following degrees were earned at Texas A&M University by individuals 
associated with this research project: 

El-banna, Hesham M., Master of Science (Aerospace Engineering), May 1988. 

El-banna, Hesham M., Doctor of Philosophy (Aerospace Engineering), August 
1992. 

Arslan, Alain, Master of Science (Aerospace Engineering), December 1993. 

The following publications resulted from research associated with this project: 

El-banna, H. M., "Numerical Computation of Aerodynamic Sensitivity 
Coefficients in the Transonic and Supersonic Regimes," Master of Science Thesis, 
Aerospace Engineering Department, Texas A&M University, College Station, Texas, 
May 1988. 

El-banna, H. M. and Carlson, L. A., "Determination of Aerodynamic Sensitivity 
Coefficients in the Transonic and Supersonic Regimes," AIAA Paper 89-0532, 
January 1989. 

El-banna, H. M. and Carlson, L. A., "Determination of Aerodynamic Sensitivity 
Coefficients Based on the Transonic Small Perturbation Formulation," Journal of 
Aircraft , Vol. 27, No. 6, June 1990, pp. 507-515. 

El-banna, H. M. and Carlson, L. A., "Determination of Aerodynamic Sensitivity 
Coefficients Based on the Three-Dimensional Full Potential Equation," AIAA Paper 
92-2671, June 1992. 

El-banna, H. M. and Carlson, L. A., "A Compendium of Transonic Aerodynamic 
Sensitivity Coefficient Data," TAMRF Rept. 5802-9203, Texas A&M Research 
Foundation, College Station, TX, July 1992. 

El-banna, H. M., "Aerodynamic Sensitivity Analysis in the Transonic Regime," 
Doctor of Philosophy Dissertation, Aerospace Engineering Department, Texas A&M 
University, College Station, Texas, August 1992. 

Carlson, L. A. and El-banna, H. M., "Determination of Aerodynamic Sensitivity 
Coefficients Based on the Three Dimensional Full Potential Equation: Users Guide 
for Analysis/Sensitivity Program and Graphics Program, " Aerospace Engineering 
Department, Texas A&M University, College Station, Texas, August 1992. 
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Carlson, L. A. and El-banna, H. M., "Determination of Aerodynamic Sensitivity 
Coefficients for Wings in Transonic Flow," Proceedings of the 3rd Pan American 
Congress of Applied Mechanics , D. T. Mook, editor, Sao Paulo, Brazil, January 
1993, pp. 13-16. 

Arslan, A. E. "Analysis and Numerical Computation of Sensitivity Derivatives in 
the Transonic Regime," Master of Science Thesis, Aerospace Engineering, Texas 
A&M University, College Station, Texas, December 1993. 

El-banna, H. M. and Carlson, L. A., "Aerodynamic Sensitivity Coefficients Using 
the 3-D Full Potential Equation," accepted for publication in the Journal Of 
Aircraft , 1994. 

Arslan, A. E. and Carlson, L. A., "Integrated Determination of Sensitivity 
Derivatives for an Aeroelastic Transonic Wing," Submitted to the 5th 
AIAA/USAF/NASA/ISSMO Symposium on Multidisciplinary Analysis and 
Optimization, September 1994. 

Copies of some of these publications are included in the appendix of this report. 
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Abstract 


The quasi-analytical approach is developed 
to compute airfoil aerodynamic sensitivity 
coefficients in the transonic and supersonic 
flight regimes. Initial investigation verifies 
the feasibility of this approach as applied to 
the transonic small perturbation residual 
expression. Results are compared to those 
obtained by the direct (finite difference) 
approach and both methods are evaluated ~to 
determine their computational accuracies and 
efficiencies. The quasi-analytical approach is 
shown to be superior and worth further 
investigation . 

Nomenclature 

Al, A2 Coordinate stretching constants 

C Maximum camber in fraction of chord 

Cp Pressure coefficient 

IM, JM Grid dimensions 
JB Row above airfoil 

Chordwise location of maximum camber 
M Mach number 

R Residual expression 

T Maximum thickness in fraction of 

chord 

XD Design variable 

f, g Cartesian coordinate stretching 

functions 

x, y Cartesian coordinates 

oc . Angle of attack 

7 Ratio of specific heats 

T Circulation 

<p Perturbation potential function 

ACp Cp x - Cp u 

Subscripts 

«=° Free stream condition 

b Body 

p Pressure 

u, 1 Upper, lower 
TE Trailing edge 

Introduction 

Over the past few years, computational fluid 
dynamics has evolved rapidly as a result of the 
immense advancements in the computational field 
and the impact of the use of computers on 
obtaining numerical solutions to complex 
problems. Accordingly, researchers are now 
capable of calculating aerodynamic forces on 
w t n g - b o dy - nacelle - empennage conf igura t ions 

subject to subsonic or transonic flows. A next 
logical step would be to compute the sensitivity 
of these forces to configuration geometry. 

* Graduate Research Asst., Student Member AIAA 
** Professor, Aerospace Engr. , Assoc. Fellow AIAA 

Copyright © American Institute of Aeronautics and 
Astronautics. Inc., 1989 All rights reserved. 


In order to Improve the design of transonic 
vehicles, design codes are being developed which 
use optimization techniques; and, in order to be 
successful, these codes require aerodynamic 
sensitivity coefficients, which are defined as 
the derivatives of the aerodynamic functions with 
respect to the design variables. Obviously, it is 
desirable that such sensitivity coefficients be 
easily obtained. Consequently, the primary 
objective of this effort is to investigate the 
feasibility of using the quasi - analyticajl. 
method^ for calculating the aerodynamic 
sensitivity derivatives in the transonic and 
supersonic flight regimes. As part of this work, 
the resulting sensitivity coefficients are 
compared to those obtained from the finite 
difference approach. Finally, both methods are 
evaluated to determine their computational 
accuracies and efficiencies. 

In the transonic regime, a variety of 
flowfield solution methods exist. These range 
from full Navier- Stokes solvers to transonic 
small perturbation equation solvers. The 
complexity of the equations that need to be 
solved depends upon the flow phenomena in 
question and the objective of the analysis. Since 
it Is not the objective of this work to develop 
flowfield algorithms, the present research uses 
the nonlinear transonic small perturbation 
equation to determine and verify efficient 
methods for calculating the aerodynamic 
sensitivity derivatives. In addition, only two 
dimensional results will be presented in this 
initial work. 

Background 

Most recently, sensitivity methodology has 
been successfully used in structural design^ and 
optimization programs^ primarily to assess the 
effects of the variation of various fundamental 
properties relative to the Important physical 
design variables. Moreover, researchers have 
developed and applied sensitivity analysis for 
analytical model improvement and assessment of 
design trends. In most cases, a predominant 
contributor to the cost and time in the 
optimization procedures is the calculation of 
derivatives. For this reason it is desirable In 
aerodynamic optimization to have efficient 
methods of determining the aerodynamic 
sensitivity coefficients and, wherever possible, 
to develop appropriate numerical methods for such 
computations . 

Currently, most methods for calculating 
transonic aerodynamic sensitivity coefficients 
are based upon the finite difference 
approximation to the derivatives. In this 
approach, a design variable is perturbed from its 
previous value, a new complete solution is 
obtained, and the differences between the new and 
the old solutions are used to obtain the 
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sensitivity coefficients. This direct, or brute 
force, technique has the disadvantage of being 
potentially very computer intensive, especially 
if the governing equations are expensive to 
solve. Accordingly, the need to eliminate these 
costly and repetitive analyses is the primary 
motivation for the development of alternative 
efficient computational methods to determine the 
aerodynamic sensitivity coefficients. 

Problem Statement 


Based on the foregoing discussion, the 
current problem is formulated starting from the 
generic quasi -analytical approach and manipulated 
according to the rules given in Appendix A of 
Ref.l for the derivation of the general 
sensitivity equation. This general sensitivity 
equation is then applied to the residual 
expression (R) of the transonic small 
perturbation equation, which Is a simple and 
adequate description of the nonlinear phenomena 
occurring in the transonic regime. Although this 
expression is nonlinear in the perturbation 
potential (<p) , the general sensitivity equation, 
Eq.(l), is linear with respect to the unknown 
sensitivity (dip/dXD^). It is to be noticed that 
the practical implementation of the above step is 
not achieved until the residual expression is 
approximated on a finite domain and the 
mathematical form of the problem rendered to that 
of one In linear algebra. This discretization 
process is explained In detail In Ref. 6. 

Thus, the quasi -analytical method, as 
applied to the residual expression of the 
transonic small perturbation equation, yields the 
sensitivity equations, 



R - (Bl + B 2 <P X ) qp xx + <Pyy - 0 
B 1 - 1 - 
B 2 — - (7+1) M*. 3 
<p ■ cp (x,y, XD) 


( 1 ) 


( 2 ) 


(3) 


and the Kutta Condition 


AP - 0 (r - A<p - const.) , x*te < x < « 

Equation (1) Is discretized into a system of 
linear equations to be solved for the unknown 
sensitivity vectors. The solution of this system 
is obtained efficiently by using either a direct 
or an iterative procedure that allows for 
multiple right hand sides. This approach is 
explained in the following section and has the 
advantage that several unknown vectors can be 
obtained simultaneously, each vector representing 
the sensitivity of the potential ( <p ) with respect 
to some design variable XD^. 

At this stage, it Is convenient to define 
the vector of design variables 


XD - { XD 1p XD 2 , - . - . XD n } (7) 

and to exactly determine which variables 
influence the solution of Eq.(2). In doing so, 
the relation between the sensitivity coefficients 
corresponding to these variables and the form of 
the optimization algorithm that utilizes this 
information needs to be considered. Notice that 
the derivatives computed in this study, namely, 
the first partial derivatives, are adequate for a 
typical optimization routine if it were to be 
applied to the present two dimensional problem. 
Notice also that some optimization studies might 
require higher derivatives. 

For the transonic flow problem, an 
appropriate choice of the first design variable 
is the free stream Mach Number (M«,) - This 
variable appears in the governing Eq.(2) 
and has an important influence on the character 
of the equation via its influence on local Mach 
number ( for M<l,the equation is elliptic, for 
M>1, the equation Is hyperbolic ) and thus on the 
nature of the solution. For this reason, it is 
desirable to have M* as one of the design 
variables . 

Next, it is appropriate to examine the 
boundary condition given by Eq.(5). In the 
transonic small perturbation formulation, the 
angle of attack (oc) enters the problem through 
the boundary condition and thus, 


r u 

1 


dy 

— - y u ' - a 

dx b 1 
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XD — set of design variables 
XD^ ■ 1^ design variable 
subject to the airfoil boundary condition, 

T d y 1 


*>y(xb,0) 


— - F(x , XD) 

dx Jb 


the infinity boundary condition, 
for M*<1 

- - r*/(2ir), 9 - nx/2, n - 0,1,2 ,3, 4 


or for M^l 


<Poo - 0 , 8 — nir/2 , n - 1 , 2 , 3 

?> x - 0 , d - nfr/2 , n - 0,4 


(A) 


( 5 ) 


For simplicity, the function (F) should be easily 
differentiable with respect to the design 
variables defining the airfoil geometry. This 
desirable feature simplifies the computation of 
the right hand side term of the sensitivity 
equation. Therefore, it would seem plausible to 
have a simple analytical expression for modeling 
the upper and lower surfaces of the airfoil. 

For the present studies, it was decided to 
limit consideration to two basic airfoil 
sections, namely parabolic -arc sections, and the 
NACA four-digit sections, whose families of wing 
sections are obtained by combining a mean line 
and a thickness distribution^. The resultant 
expressions possess the necessary features that 
suit the problem, mainly the concise description 
of the airfoil surfaces in terms of several 
geometric design variables. The expressions are 
as follows : 




For parabolic -arc sections 
2 v 2 


fC(2Lx— x 2 )/L 

yuH 


±2Tx(l-x) , x±;L 


flci (l-L)+2Lx-x 2 ]/(l-L) 2 ±2Tx(l-x) , x>L 


(9) 


For NACA four-digit sections 

fC(2Lx-x 2 )/L 2 ± 5T(0.29697x-0. 126x 
I -0 . 3516x+0 . 2843x -0.1015x ), x<L 


y ;1 


C[ (1-2L) +2Lx-x ? ]/ (1— L) 2 + 5T(0.29697x 
-0 . 126x-0 . 3516x +0 . 2843x -0.10l5x ), x>L 


( 10 ) 


Each of the quantaties C, L, and T is 
expressed as a fraction of the chord (e.g. if T 
is 6% chord then T - 0.06). Differentiating 

Eqs . (9) and (10) with respect to x and 
substituting the results into Eq.(8) yields : 


For parabolic-arc sections 

F u ( i - 2C(L-x)/LL ± 2T(l-2x) - « (11) 

For NACA four-digit sections 

F u<1 - 2C(L-x)/LL ± 5T(0.148A5/7x-0.126 

-0.7032x+0.8529x 2 -0.406x S ) - <* (12) 


where 


L 2 , x^L 
(1— L) 2 , x>L 


(13) 


f - <d{/dx) - A 2 (l-« 2 ) 
g - (d,/dy) - Aid-, ) 

so that, 

<p x - f<P£ 

<fy " g <Pn 
<Pxx " f ( f ^>^ 

<Pyy ~ B 

Substituting from Eqs . (19) - (22) into Eq.(2), 
yields the transformed residual expression. 

R - [B 1 +B 2 f^] f(f<P{)f + g(g^)r? - 0 ( 23 > 

This equation is solved numerically by an 
approximate factorization scheme 9 in which the 
objective is to force the residual to zero at 
each point of the computational domain. In finite 
difference form, Eq.(23) can be written as, 

R i,j * l B l + ^(Vi+l , j-^i-1 , j )/( 2a O 1 

[‘'i, j f i+b('«Pi+l, j”^i, j) 

-<2i/£ # j-1) f t j -qPi-l , j ) 

, j ) f i-3/2^i-l , j^i-2 . j > 1 
+ [gj+H<<Pi,j+r^i,j> 

2 

-gj_lj(Vi, j^i,j-l>] gj/ A 9 ( 24 > 

where 


(17) 

(18) 


(19) 

( 20 ) 
( 21 ) 
(22) 


vi j - 1 if point (i.j) is subsonic 

- 0 if point (i, j ) is supersonic 


Eqs. (11) and (12) are simple analytical 
expressions in terms of the four variables T, L, 
C, and a. Thus, 

XD - { T, Meo, «, L f C } (14) 

represents the complete set of design variables 
that define the present two-dimensional airfoil 
sensitivity problem. Notice that these variables 
are completely uncoupled and, thus the 
sensitivity equation can be solved independently 
with respect to each variable®. 

Mathematical Treatment and Solution Procedure 

Problem Discretization 

Equation (1) represents the general 
sensitivity equation applied to the residual R. 
Now, in order to solve the problem numerically, 
Eq.(2) is formulated computationally on a finite 
domain. This transformation is achieved by using 
a stretched Cartesian grid that maps the infinite 
physical domain onto a finite computational grid. 
In this study, the grid used is based upon a 
hyperbolic tangent transformation that places the 
outer boundaries at infinity. Accordingly, the 
computational variables used are given by, 

£ - tanh A 2 x (15) 

rj - tanh A^y ( 16 ) 

In addition, the stretching functions are 
defined as. 


Eq.(24) is the discretized form of the residual 
at a general point (i.j) in terms of <p values at 
surrounding points. Consequently, R at i,j can be 
viewed as a function of the <p values at 
neighboring points; and, therefore, the 
differentiation of the residual expression is 
straight forward. 

Differentiation of the Residual 


Rearranging Eq.(24) yields 

R i,j “ c m,j + c 2^i+l,j < Pi-l,j + c 3*l+l,J*l,J 

+ c 4?i-l, j^i, j + c 5*i+l, j<Pi-2, j 

2 2 

+ c 6Pi-l,j<Pi-2,j + c 7n-l,j + c 8<Pi+l,j 

+ c 9 < Pi+l,j + c 10 < Pi-l,j + c im,j+l 

+ c 12 <Pi,j-l + c 13<Pi-2,j (25) 

For a fixed computational grid, the 

coefficients c^, c 2 , ... , c ^3 are functions only 

of B]_ and B 2 which in turn are functions of M^. 
This fact is used when differentiating Eq.(25) 
with respect to in order to obtain the right 
hand side (aR/aM*,) . It is also necessary to 
consider the treatment of various types of grid 
points and examine the effect on the general 
residual expression. Several groups of points, 
such as those adjacent to the airfoil, to the 
wake cut, and to infinity boundaries, need 
special treatment. Accordingly, it is necessary 


if 
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to revise the residual expression at these 
boundary points to include the boundary 
conditions. The resulting updates are used to 
modify the residual equation, Eq.(25), and yield 
a set of expressions, each being valid for a 
group of boundary points. The details of these 
operations are found in Ref. 6. 

In setting up the complete quasi -analytical 
problem the circulation and its dependence upon 
trailing edge potentials must be carefully 
included. Since the circulation is determined by 
the difference in potentials at the trailing 
edge , 


r "" <?uTE “ SPITE 


(26) 


or, by interpolating the trailing edge values 


T - T X ( 1.5 

- 0.5 

* ( X.5 

- 0.5 


tolTE-1, JB ”^ITE-l t JB-l) 
^ITE-l,JBfl^PlTE-l,JB-2) 1 
(<PITE,JB "^ITE.JB-l > 
(^ITE.JB+1 ^ I TE , JB- 2 ) 1 


where 


(27) 


T 2 - [ f(x-0.5) - £(ITE-1) ] / Af (28) 

Ti - [ 1 - T 2 ] (29) 


and since a branch cut extends from the trailing 
edge to downstream infinity, the trailing edge 
potentials appear in the residual expressions for 
points along the branch cut. In addition, since 
in the two dimensional case the infinity boundary 
conditions are proportional to the circulation, 
the trailing edge potentials also appear in the 
residual expressions at points adjacent to the 
outer boundaries. Consequently, the- resultant 
matrix (3R /d<p) , while banded, also contains many 
nonzero elements far from the central band. 
Notice that the presence of these elements 
greatly complicates the rapid and efficient 
solution of the sensitivity equation, Eq.(l). 

The resulting residual expressions are 
differentiated analytically with respect to the 
potential (<p) . To be more specific, each equation 
is differentiated with respect to the potential 
at neighboring points and trailing edge points 
(the later enters as a result of the implicit 
nature of the circulation effects). These points 
are denoted by the counters (ii.jj) and are given 

b y . 


(i.j-1). (i.j), (i,j+D, (i-2 , j ) , (1-1J), 

( i+1 , j ) , (ITE-1 , JB-2) , (ITE-l.JB-l), (ITE- 

1,JB), (ITE-1, JB+1) , (ITE.JB-2), (ITE.JB-l), 

( ITE , JB) , (ITE.JB+l). 

Solution about a Fixed Design Point 

Once the residual relations are obtained, 
the actual coefficients are assembled by 
evaluating the appropriate analytical 
expresssions using a flowfield solution obtained 
from Eq.(2) for a given set of conditions (i.e. 
about a fixed design point). Similarly, the right 
hand sides are evaluated by differentiating the 
analytical expressions for the residual with 
respect to each design variable. Again, the 
details and results of these steps are found in 


Ref. 7. 

The end result is that the coefficient 
matrix (3R.£ j/<^Pii jj) is of size (IM— 2)*(JM— 
2)x(IM-2)*(jif-2) for a general (IM*JM) grid. This 
system is large, of block structure, diagonally 
dominant, and sparse; and, while banded, also 
contains many nonzero elements far from the 
central band. As a result of this size and 
structure, it is obvious that a reasonably fast 
scheme for solving Eq.(l) is needed. 

Currently, it is very difficult to single 
out an optimum routine that handles a general 
large sparse system of linear equations for which 
the coefficient matrix is unsymmetric. This is 
due to the fact that, unlike the theory of 
symmetric matrices, the theory of general 
unsymmetric matrices is more involved and has yet 
to be developed. Since research in the above 
areas is currently very active and specialized, 
any attempt to cover these topics in detail would 
be laborious. For this reason, it was decided to 
use a few general approaches that were available 
in the literature and that could be integrated 
into the sensitivity codes with adjustments. This 
approach would allow an evaluation of the overall 
cost involved in solving the current two- 
dimensional problem and would give a crude 
estimate of the effort Involved in solving a 
three dimensional problem. 

The first solver is based on standard 
Gaussian Elimination with partial pivoting and 
full storage. The second is based on triangular 
decomposition 1 -® and uses a compact storage scheme 
that avoides handling the zero entries and 
therefore should be more efficient than standard 
Guasslan Elimination. The third solver is based 
on a Gauss -Seidel iterative scheme 11 and was not 
optimized for speed (through the choice of 
optimum acceleration parameters) but uses sparse 
matrix technology in processing only the nonzero 
elements. The fourth and last solver used is 
based on the conjugate gradient method 12 . 
Handling the sparsity pattern for the third and 
fourth solvers is achieved by assembling the 
symbolic part of the coefficient matrix only once 
for a given grid size and given free -stream 
(subsonic versus supersonic). The resultant 
structure is then stored on a diskfile. Before 
the numerical part is executed, the symbolic 
information is read into the code and used 
directly to assemble the new matrix. This 
procedure is followed in order to reduce the time 
consumed in assembling the coefficient matrix. 
Notice also that in the Gauss-Seidel and 
conj ugate - grad ient solvers that the error 
tolerances for the coefficients involving maximum 
thickness, free stream Mach number, and location 
of maximum camber were l.E-06 while those on 
angle of attack and maximum camber were 1 . E-04 . 

Once the sensitivities of the potentials, 
and thus the Cp distribution, to the design 
variables are known, the sensitivity of the lift 
coefficients to the design variables can be 
easily computed. To minimize errors, these 
coefficients are computed using 

C L - 2 r - 2 (<P u te~^ 1TE) (30) 
and hence , 

3Cl/ 3XD^ - 2 (3<p u j£/aXD^-J<p^ TE /aXD^) (31) 



Finally, all methods used for computing the 
derivatives are compared to the finite - 
difference approach and the results are presented 
and evaluated to determine the computational 
accuracy and efficiency (with regards to time) of 
each method. 

Test Cases 

In this study, the quasi -analytical method 
has been used to determine the aerodynamic 
sensitivity coefficients at three freestream 
Mach numbers (M^ - 0.2, 0.8, 1.2) for two 
arbitrarily selected airfoils, each at one degree 
angle of attack. The first is a cambered 
parabolic arc section having 1% camber at 40% 
chord, a maximum thickness of 6% at 50% chord, 
and which is designated P1406; and the second is 
a NACA 1406 airfoil. Since most of the 
interesting captured phenomena were found to be 
identical for both airfoils, only results for the 
NACA 1406 airfoil are presented in this paper. 

In the following, two types of results will 
be presented. The first will be plots of Cp 
versus chord for the three chosen Mach numbers. 
The second will be the corresponding plots of 

(3c p/3T), (acp/ahu), (acp/aac), (acp/ac), and 

(aep/aL) obtained by the quasi -analytical method. 
In addition, all of the figures will also contain 
results obtained using the direct (finite 
difference) approach in which each design 
variable was individually perturbed by a small 
amount, typically 0.001, and a new flowfield 
solution obtained. Then the sensitivities were 
computed using ACp/AXD and are shown via dashed 
lines. In many cases the lines are coincident 
with the quas i- analytical results and cannot be 
observed. Table I compares results obtained by 
the two methods, and in most cases the agreement 
is within one percent. 

In all cases, an 81*20 stretched Cartesian 
grid was utilized. In addition, for these 
studies, the flowfield was normally computed 
using double precision arithmetic and the maximum 
residual reduced eight orders of magnitude. It 
was felt that this level of convergence was 
necessary in order to accurately evaluate 
sensitivity coefficients using a finite 
difference approach, although such convergence 
may not be required in the flowfield solver for 
the quasi - analytical method. 

Results and Discussion 
Subsonic Case (M^ - 0.2) 

Initial studies concentrated on subsonic 
cases since at least approximate results would be 
known from thin airfoil theory. Figure 1 shows 
the pressure distribution for the NACA 1406 
airfoil while Figs. 2a and 2b show the sensitivity 
of the pressure to thickness for the same 
airfoil. As expected from thin airfoil theory, 
the upper and lower surface values are 
essentially identical and the difference is very 
small everywhere. Also shown on the same figure 
(and on subsequent figures) by the dashed line is 
the result obtained by using the finite 
difference approach; and as can be seen, the 
agreement betweeen the two approaches is 
excellent. 

The sensitivity of pressure to freestream 


Mach number is plotted on Figs. 3a and 3b. It is 
noticed that while the profiles for the upper and 
lower surfaces are similar, they are not equal in 
magnitude, indicating a nonlinear variation with 
Mach number as predicted by simple Prandtl- 
Glauret Theory. However, as indicated by the 
results on Fig. 3b, the magnitudes for this 
subsonic Mach number are very low. 

The sensitivity of the pressure coefficients 
to angle of attack are shown for this case on 
Figs. 4a and 4b. As expected from linear thin 
airfoil theory, the upper and lower surface 
curves are essentially equal in magnitude but of 
opposite sign. Not surprisingly, the sensitivity 
of the delta Cp variation, Fig. 4b, has the shape 
of the pressure difference curve for a flat plate 
at angle of attack; and its magnitude, 
particularly near the leading edge is quite 
large . 

On Figs. 5a and 5b is plotted the sensitivity 
of the pressure coefficient to the amount of 
maximum camber. Since camber contributes to lift, 
it is expected from thin airfoil theory that 
these values should be "equal but opposite in 
sign" for the upper and lower surfaces. In 
addition, the pressure difference curve has the 
correct shape for that associated with a 14 mean 
line with the peak occur ing at 30% chord 7 and has 
magnitude comparable to those for the (3Cp/3«) 
curves . 

Finally, the sensitivity of pressure to the 
location of the maximum camber point is portrayed 
on Figs. 6a and 6b, and to say the least the 
results are interesting. Since maximum camber 
location affects the camber profile and hence 
lift, the equal and opposite behavior of the 
upper and lower surface coefficients is expected. 
In addition, the pressure difference sensitivity 
is primarily negative forward of the point of 
maximum camber and positive aft of it. This 
result indicates that if the location of maximum 
camber were moved rearward slightly (i.e. a 
positive AL) that lift would be decreased on the 
forward portion of the airfoil and increased on 
the aft portion of the airfoil, which is in 
agreement with the results presented in Ref. 7. 

Transonic Case (M^ - 0,8) 

At — 0.8, the flow about the NACA 1406 
airfoil has a strong shock at 40% chord, Fig. 7; 
and the lower surface is entirely subcritical. As 
a consequence, the variation with chord of the 
sensitivity coefficients is considerably 
different than in the subsonic case. 

Figs. 8a and 8b show the sensitivity of 
pressure to the maximum thickness; and while the 
lower surface profile is similar to that obtained 
at subsonic conditions, the upper surface curve 
and the pressure difference coefficient plot show 
the effect of the upper surface shock wave. The 
large peak on the curves corresponds to the 
location of the shock wave and indicates that the 
shock wave location is very sensitive to maximum 
thickness. Notice on Figs. 8a and 8b the excellent 
agreement of the quasi -analytical results 
Indicated by the solid lines with those obtained 
using the finite -differece approach (dashed 
lines) . 

The results for (3Cp/3M^), which are shown 
on Figs. 9a and 9b, are similar. The lower surface 
curve is typical of a subsonic flow, while the 
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upper surface and the pressure difference 
coefficients reflect the presence of the upper 
surface shock wave. Similar comments can be made 
for the remaining design variable coeficients, 
which are plotted on Figs. 10, 11, and 12. 

Examination of the curves in the vicinity of 
the shock wave location indicates that the 
pressure sensitivity and indirectly the shock 
wave location is about equally influenced by the 
maximum thickness, freestream Mach number, and 
angle of attack. However, in comparison it is 
relatively insensitive to location of maximum 
camber; but, perhaps surprisingly so, the 
pressure is twice as sensitive to the amount of 
maximum camber as it is to the other design 
variables. It should also be noticed that the 
lift Is most sensitive to angle of attack and to 
maximum camber. 

In addition, Fig. 11 shows a discrepancy 
between the results obtained by the direct 
approach and those obtained thru the quasi - 
analytical method. It will be shown in the 
following section that this discrepancy is 
related to the choice of the step size used in 
computing the finite -difference solution, thus 
revealing a significant deficiency in computing 
the sensitivity derivatives in nonlinear regimes 
via the finite-difference approach. 

Supersonic Case - 1.2) 

In order to investigate the applicability of 
the quasi - analytical method at supersonic 
freestream Mach numbers, solutions were obtained 
for the NACA 1406 airfoil at Mach 1.2. At this 
condition, the flow is transonic in that the bow 
shock is detached, and there is a region of 
subsonic flow extending to approximately the 
quarter chord. Fig. 13. Figures 14-18 show the 
pressure sensitivities for these cases, and Table 
I lists the lift sensitivities. 

Examination of the plots shows that the 
pressure sensitivity coefficients have different 
trends and magnitudes from those computed for 
subsonic freestream supercritical conditions, and 
that they are approaching the form expected from 
supersonic linear theory. These changes are 
particularly evident in the lift derivatives 
presented in Table I. Notice that the derivatives 
with respect to the design variables maximum 
thickness, Mach number, and location of maximum 
camber have switched sign. In addition, as 
expected from linear theory, the influence of 
camber on lift has decreased significantly; and 
at Ken - 1.2 is only about 13% of the angle of 
attack effect as compared to a factor of about 
two at Moj - 0,8. Notice also that Fig. 15 shows a 
discrepancy similar to that found in Fig. 11. 

Time Comparisons 

Obviously, in the development of the quasi - 
analytical method It was hoped that not only 
would this approach yield accurate values for the 
aerodynamic sensitivity coefficients but also 
that It would be more efficient than the brute 
force finite difference approach. Table II 
presents some comparisons concerning the amount 
of computational effort required to obtain 
solutions by the two approaches. 

In comparing the values, several items 
should be kept in mind. First, It has been 


assumed that the finite difference approach will 
require six independent solutions. In practice it 
might be possible to start each finite difference 
solution from a previous solution and, thus, 
decrease the time to convergence. However, to be 
accurate, the finite difference approach will 
probably require double precision and will have 
to be extremely well converged ( i . e . 1 . E-08) . 
Nevertheless, the values for the finite 
difference approach probably should be viewed as 
maximum values. 

Second, the methods used for obtaining the 
sensitivity coefficients have not been optimized 
and, as mentioned earlier, may not even be 
optimum; and the flowfield solution required for 
the quasi -analytical approach may not need double 
precision and may not have to be as tightly 
converged. Thus the values shown for the quasi - 
analytical approach should also be viewed as 
maximum values. 

In spite of these limitations, results 
obtained by direct methods do indicate, that the 
quas i - analytical method is at transonic 
conditions potentially more computationally 
efficient than the brute force finite difference 
approach . 

Notice that in this study, the initial guess 
used in computing the sensitivity derivatives via 
iterative methods was arbitrarily chosen as the 
zero-vector. In addition, time comparisons 
presented in Table II show that iterative methods 
are In general less efficient than direct methods 
if the derivatives for the current two- 
dimensional problem were sought about some 
general design point. However, if the objective 
Is to incorporate the sensitivity derivatives In 
an optimization loop (i.e. to use the derivatives 
in a continuation problem), then, a good initial 
guess (which in this case would be available) 
would enhance convergence and the overall cost of 
computing the derivatives using iterative methods 
might be reduced. These points should be taken 
into consideration when a sensitivity study is to 
be integrated into an optimization procedure. 

Additional Test Cases 

The first group of cases are carried out to 
investigate the performance of the NACA 1406 
airfoil for a range of Mach numbers from 0.79 to 
0.86 in increments of 0.01. As shown in Fig. 19, 
this range of transonic Mach numbers encompasses 
the development of the shock wave on the upper 
surface of the airfoil. Also, as shown on Figs. 20 
and 21 for the cases involving thickness, Mach 
number, and maximum camber, the quasi -analytical 
derivatives are In the vicinity of the shockwave 
frequently different from those obtained by the 
f inite -difference approach. This discrepancy 
raises two questions -- Vhat is the cause of the 
disagreement and which set of derivatives is more 
accurate ? Examination of the variation of the 
integrated coefficient, dCL/dXD^ with M^, which 
is portrayed on Fig. 22, shows that the quasi - 
analytical results are smooth and follow a 
definite trend while the finite difference values 
are at best "discontinuous". Consequently, it is 
concluded that the f lnlte - difference results are 
less accurate. In order to observe the 
performance of the finite -difference approach in 
the transonic regime, it is necessary to examine 
the effect of changing the step size (delta of 


the design variable) on the computed derivatives. 
Four different values for the step size (l.E-03, 

l.E-04, l.E-05, and l.E-06) were chosen and 
applied to the NACA 1406 at a Mach number of 
0.84. Examination of the results (Table III) show 
that as the step size is decreased, the finite 
difference lift coefficient sensitivity 
derivatives approach the values computed by the 
quas i - analytical method. However, In some cases, 
for small AXDj^ values, oscillations in the 
pressure coefficient sensitivity derivatives have 
been observed depending upon the machine used and 
the method of storing and retrieving the data. 
These oscillations combined with the difficulty 
of properly choosing a suitable finite difference 
AXD^ a priori indicates that the finite 
difference approach is probably not a practical 
method of efficiently computing sensitivity 
coefficients. On the other hand, the present 
results demonstrate that the quas i -analytical 
method can be used accurately to obtain such 
coefficients in the transonic flight regime. 

Conclusion 

Based upon these investigations and results, 
it is concluded that the quasi- analytical method 
Is a feasible approach for accurately obtaining 
transonic aerodynamic sensitivity coefficients in 
two dimensions. The results obtained from the 
quasi -analytical method are almost Identical to 
those obtained by the brute force (finite 
difference) technique. Furthermore, the study 
indicates that obtaining the quasi -analytical 
transonic derlvates using a direct solver is more 
efficient than computing the derivatives by the 
finite difference method. 
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Determination of Aerodynamic Sensitivity Coefficients 
Based on the Transonic Small Perturbation Formulation 
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The quasianatytical approach is developed lo compute airfoil aerodynamic sensitivity coefficients in the 
transonic and supersonic night regimes. Initial investigation verifies the feasibility of this approach as applied 
to the transonic small perturbation residual expression. Results are compared to those obtained by the direct 
(finite difference) approach, and both methods are evaluated to determine their computational accuracies and 
cffic.enc.es. The quasianalytieal approach is shown to yield more accurate coefficients and is potentially more 
efficient and worth further investigation. 
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Nomenclature 

Id 

C 

= maximum camber in fraction of chord 

M 

Cp 

= pressure coefficient 


IMJM 

= grid dimensions 


JB 

= row above airfoil 

! = 

L 

= chordwise location of maximum camber 

— 

M 

= Mach number 


R 

= residual expression 


T 

= maximum thickness in fraction of chord 


XD 

= design variable 

_ 

/.* 

= Cartesian coordinate stretching functions 



= Cartesian coordinates 



= computational variables 


a 

= angle of attack 


7 

= ratio of specific heats 


r 

= circulation 


<p 

= perturbation potential function 


A Cp 

= Cp x - Cp u 


Subscripts 


oo 

= freestream condition 

„ 

b 

— body 

L- 

P 

= pressure 


u t 1 

= upper, lower 


TE 

= trailing edge 


Introduction 

O VER the past few years, computational fluid dynamics 
has evolved rapidly as a result of the immense advance- 
ments in the computational field and the impact of the use of 
computers on obtaining numerical solutions to complex prob- 
lems. Accordingly, researchers are now capable of calculating 
aerodynamic forces on wing-body-nacelle-empennage config- 
urations. A next logical step would be to compute the sensitiv- 
ity of these forces to configuration geometry. 

In order to improve the design of transonic vehicles, design 
codes are being developed that use optimization techniques; 
and, in order to be successful, these codes require aerodynamic 
sensitivity coefficients, which are defined as the derivatives of 
the aerodynamic functions with respect to the design variables. 
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Obviously, it is desirable that such sensitivity coefficients be 
easily obtained. Consequently, the primary objective of this 
effort is lo investigate the feasibility- of using the quasi- 
analytical method 1 3 for calculating the aerodynamic sensitiv- 
ity derivatives in the transonic and supersonic flight regimes. 
As part of this work, the resulting sensitivity coefficients are 
compared to those obtained from the finite difference ap- 
proach. Finally, both methods are evaluated to determine their 
computational accuracies and efficiencies. 

In the transonic regime, a variety of flowfield solution meth- 
ods exist. These range from full Navier-Stokes solvers to tran- 
sonic small perturbation equation solvers. The complexity of 
the equations that need to be solved depends upon the flow 
phenomena in question and the objective of the analysis. Since 
U is not the objective of this work to develop flowfield algo- 
rithms, the present research uses the nonlinear transonic small 
perturbation equation to determine and verify efficient meth- 
ods for calculating the aerodynamic sensitivity derivatives. In 
addition, only two-dimensional results will be presented in this 
initial work. 


Background 

Most recently, sensitivity methodology has beensuccessfulty 
used in structural design 2 and optimization programs 3 primar- 
ily to assess the effects of the variation of various fundamental 
properties relative to the important physical design variables. 
Moreover, researchers have developed and applied sensitivity 
analysis for analytical model improvement and assessment of 
design trends. In most cases, a predominant contributor to the 
cost and time in the optimization procedures is the calcula- 
tion of derivatives. For this reason, it is desirable in aero- 
dynamic optimization to have efficient methods of determin- 
ing the aerodynamic sensitivity coefficients and, wherever pos- 
1 sible, to develop appropriate numerical methods for such 
computations. 

Currently, most methods for calculating transonic aero- 
dynamic sensitivity coefficients are based upon the finite dif- 
ference approximation to the derivatives. In this approach, a 
design variable is perturbed from its previous value, a new 
complete solution is obtained, and the differences between the 
new and the old solutions arc used to obtain the sensitivity 
coefficients. This direct, or brute force, technique has the 
disadvantage of being potentially very computer intensive, es- 
pecially if the governing equations are expensive to solve. In 
addition, it is difficult to guarantee the accuracy of (he deriva- 
tives obtained by the finite difference method. Accordingly, 
the need to eliminate these costly and repetitive analyses is the 
primary motivation for the development of alternative, effi 
cicnt computational methods to determine the aerodynamic 
sensitivity coefficients. 
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Problem Statement 

. , a Based on th<£ foregoing discussiorf; the current problem is 

* 'formulated starting from the generic quasianalytical approach 
and manipulated according to the rules given in Appendix A of 
Ref. I for the derivation of the general sensitivity equation. 
This general sensitivity equation is then applied to the residual 
expression R of the transonic small perturbation equation, 
which is a simple and adequate description of the nonlinear 
phenomena occurring in the transonic regime. Although this 
expression is nonlinear in the perturbation potential <p , the 
general sensitivity equation, Eq. (1), is linear with respect to 
the unknown sensitivity (d<p/dXD;). It is to be noticed that the 
practical implementation of the above step is not achieved 
until the residual expression is approximated on a finite do- 
main and the mathematical form of the problem rendered to 
that of one in linear algebra. This discretization process is 
explained in detail in Ref. 4. 

Thus, the quasianalytical method, as applied to the residual 
expression of the transonic small perturbation equation, yields 
the sensitivity equations. 


M [ dv | _ 

d R ] 


Uw>J 


where 

R S (Z?j + B 2 <Px)<Pxx + <Pyy = 0 ( 2 ) 

B, = l - Ml 

... B 2 = ~{y+\)Ml 

<p <p(x t y t XD) ( 3 ) 

XD set of design variables 
XDi 5= ith design variable 
subject to the airfoil boundary condition, 

vv(*6» 0) = 

the infinity boundary condition, for 

V?oo — — r^/(27r), 8 = mr/2 t n =0,1, 2, 3,4 

or for > 1 

iPo, = 0, 0 = n 7t/2, n — 1,2,3 

^ = 0, 6 = n tf/2, n = 0,4 (5) 

and the Kutta condition 

AP = 0 (T = Aip = const), x T£ <x<oo (6) 


dx 


= F(x y XD) 


(4) 


and to exactly determine which variables influence the solution 
of Eq. (2). In doing so, the relation between the sensitivity 
coefficients corresponding to these variables and the form of 
the optimization algorithm that utilizes this information needs 
to be considered. Notice that the derivatives computed in this 
study, namely, the first partial derivatives, arc adequate for a 
typical optimization routine if it were to be applied to the 
present two-dimensional problem. Notice also that some opti- 
mization studies might require higher derivatives. 

For the transonic flow problem, an appropriate choice of 
the first design variable is the freestream Mach number (Af«)- 
This variable appears in the governing Eq. (2) and has an 
important influence on the character of the equation via Its 
influence on local Mach number (for M < I, the equation is 
elliptic, for Af > 1, the equation is hyperbolic) and thus on the 
nature of the solution. For this reason, it is desirable to have 
as one of the design variables. 

Next, it is appropriate to examine the boundary condition 
given by Eq. (4). In the transonic small perturbation formula- 
tion, the angle oT attack (a) enters the problem through the 
boundary condition and thus. 



For simplicity, the function F should be easily differentiable 
with respect to the design variables defining the airfoil geom- 
etry. This desirable feature simplifies the computation of the 
right side term of the sensitivity equation. Therefore, it would 
seem plausible to have a simple analytical expression for mod- 
eling the upper and lower surfaces of the airfoil. 

For the present studies, it was decided to limit consideration 
to one basic airfoil section, namely the NACA four-digit sec- 
tion, whose families of wing sections are obtained by com- 
bining a mean line and a thickness distribution. 5 The resultant 
expressions possess the necessary features that suit the prob- 
lem, mainly the concise description of the airfoil surfaces in 
terms of several geometric design variables. The expressions 
are as follows: 


c „ 

C{2Lx — x 2 )/ L 2 ± 5 T (0.2969\/x —0.1 26.x 

— 0.35 16.x 2 + 0.2843.x 3 — 0.101 5x 4 ), x < L 




* C[(\-2L) + 1Lx-x 2 ]/(\-L) 2 

± 5 7'(0.2969'/x - 0. 1 26x - 0.35 16x 2 
+ 0.2843.X 3 — 0. 101 5 jr'*). x>L 


(9) 


Each of (he quantities C, L, and T is expressed as a fraction 
of the chord. Differentiating Eq. (9) with respect to x and 
substituting the result into Eq. (8) yields 

F uA = 2 C(L -x)/LL ± 57'(0.14845/v / x —0.126 

— 0.7032.x + 0.8529.x 2 —0.406jr 3 ) - a (10) 

where 


Equation (l) is discretized into a system of linear equations 
to be solved for the unknown sensitivity vectors. In carrying 
out this step, the expressions for both the right side vector and 
the left side matrix are generated analytically. The solution of 
this system is obtained efficiently by using either a direct or an 
iterative procedure that allows for multiple right sides. Th is 
approach is explained in the following section and has the 
advantage that several unknown vectors can be obtained si- 
multaneously, each vector representing the sensitivity of the 
potential v> with respect to some design variable X D t , 

At this stage, it is convenient to define the vector of design 
variables 


L 2 , x < L 
(1 -L)\ x>L 


( 11 ) 


Eq. (10) is a simple analytical expression in terms of the four 
variables T, L , C, and a. Thus, 


XD = (7\ a, L , Cj (12) 

represents the complete set of design variables that define the 
present two-dimensional airfoil sensitivity problem. Notice 
that these variables arc completely uncoupled; and, thus, the 
sensitivity equation can be solved independently with respect 
to each variable/’ 


XD [A7}|,A7} 2i . . ,A7>„] 
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Mathematical Treatment and Solution Procedure 

Problem Discretization 

Equation (I) represents the general sensitivity equation ap- 
plied to the residual R. Now, in order to solve the problem 
numerically, Eq. (2) is formulated computationally on a finite 
domain. This transformation is achieved by using a stretched 
Cartesian grid that maps the infinite physical domain onto a 
finite computational grid. In this study, the grid used is based 
upon a hyperbolic tangent transformation that places the outer 
boundaries at infinity. Accordingly, the transformed residual 
expression is given by 

K 3 [5, + + g(g = 0 (13) 

This equation is solved numerically by an approximate factor- 
ization scheme. 7 In finite-difference form, Eq. (13) can be 
written as 

R,j ^ [fi 1 + a 2 (y> l > 1J -^- U /(2A{)]y;/A£ 2 

^ \ v iJ ft ♦ Vi (v 7 / + 1 J ~~ *Pi,j) 1 )fi — Vi {tPiJ *Pi—\ J ) 

“ (1 — ~ •Pi-lj)] 

+ [zj *■ Yi{<Pij + 1 — <Pi,j) - Zj - VlifPij - <Pi.j - »)] Zj /Ai? 2 ( 1 4) 

where 

Vi j =1 if point (ij) is subsonic 

Vi j = 0 if point (ij) is supersonic 

Eq. (14) is the discretized form of the residual at a general 
point (f\y) in terms of <p values at surrounding points. Conse- 
quently, R at / J can be viewed as a function of the <p values at 
neighboring points; and, therefore, the differentiation of the 
residual expression is straightforward. 

Differentiation of the Residual 

Rearranging Eq. (14) yields 

R t.j ~ c \<PiJ + c l<Pi + IJ + CyPi+\J<Pi.j + c 4<P,-\,j<Pij 

+ c 5<Pi+ + C#pi-\J<Pi-2J + c T<Pi-\j 2 + C 8SP,>|,; 2 

+ C*Pi+lJ + CltfPi-lJ + CuViJ+X + c \VPiJ-\ + Cn<Pi-2J 

(15) 

The coefficients Ci,c 2f . . . ,c, 3 are functions only of the 
stretching factors and of 5, and B 2 , which are functions of 
M w . This fact is used when differentiating Eq. (15) with re- 
spect to in order to obtain the right side (dR /dMJ). It is 
also necessary to consider the treatment of various types of 
grid points and examine the effect on the general residual 
expression. Several groups of points, such as those adjacent to 
the airfoil, to the wake cut, and to infinity boundaries, need 
special treatment. Accordingly, it is necessary to revise the 
residual expression at these boundary points to include the 
boundary conditions. The resulting updates are then used to 
modify the residual equation, Eq. (15), and to yield a set of 
expressions, each being valid for a group of boundary points. 
The details of these operations and the expressions for the 
coefficients C|-Cn are found in Ref. 4. 

In setting up the complete quasianaiytical problem, the cir- 
culation and its dependence upon [railing-edge potentials must 
be carefully included. Since the circulation is determined by 
the difference in potentials at the trailing edge. 


or, by interpolating the trailing-edge values 
T = Tt[l.5(<p| T n_i ys —*P\t£-\jb~ i) 

— 0.5(v?jTn- i./fl ♦ i “‘V’ire- i./a- 2 )] 

+ T 2 [ 1 .5(y?,-n- jb — jB~\i 

~ 0-5(y> ITC J8+ i — ^>, TE jb-i)\ (17) 

where 

T 2 = [a»r=0.5)-£(ITE-l)]/A£ (18) 

7-, = [1-r 2 l (19) 

and since a branch cut extends from the trailing edge to down- 
stream infinity, and trailing-edge potentials appear in the 
residual expressions for points along the branch cut. In addi- 
tion, since in the two-dimensional case the infinity boundary 
conditions are proportional to the circulation, the trailing-edge 
potentials also appear in the residual expressions at points 
adjacent to the outer boundaries. Consequently, the resultant 
matrix (3 R/dtp), while banded, also contains many nonzero 
elements far from the central band. Notice that the presence of 
these elements greatly complicates the rapid and efficient solu- 
tion of the sensitivity equation, Eq. (1). 

The resulting residual expressions are differentiated analyti- 
cally with respect to the potential <p. Specifically, each equa- 
tion is differentiated with respect to the potential at neighbor- 
ing points and trailing-edge points. The latter enter as a result 
of the implicit nature of the circulation effects. These points 
are denoted by the counters (iijj) and are given by 

O'J-1). ('.j). (ij+ I), (i-2J), (/ — l,y), (/ + IJ) 
(ITE- 1,75-2), (ITE- 1,75-1), (ITE-1,75) 
(ITE- 1,75 + 1), (ITE.75-2), (ITE, 75 - 1) 
(ITE,75), (ITE, 75 + 1) 


Solution about a Fixed Design Point 

Once the residual relations are obtained, the actual coeffi- 
cients are assembled by evaluating the appropriate analytical 
expressions using a flowfield solution obtained from Eq. (2) 
for a given set of conditions (i.e., about a fixed design point). 
Similarly, the right sides are evaluated by differentiating the 
analytical expressions for the residual with respect to each 
design variable. Again, the details and results of these steps are 
found in Ref. 4. 

The end result is that the coefficient matrix (35, j/d<Pn jj) is 
of size (IM - 2) x (JM - 2) x (JM - 2) x (JM - 2) for a general 
(JM x JM) grid. This system is large, of block structure, diag- 
onally dominant, and sparse and, while banded, also contains 
many nonzero elements far from the central band. As a result 
of this size and structure, it is obvious that a reasonably fast 
scheme for solving Eq. (1) is needed. 

Currently, it is very difficult to single out an optimum rou- 
tine that handies a general, large, sparse system of linear equa- 
tions for which the coefficient matrix is unsymmetric. This is 
because, unlike the theory of symmetric matrices, the theory of 
general unsymmetric matrices is more involved and has yet to 
be developed. Since research in the above areas is currently 
very active and specialized, any attempt to cover these topics 
in detail would be laborious. For this reason, it was decided to 
use a few general but not necessarily the most efficient ap 
proaches that were available in the literature and that could be 
integrated into the sensitivity codes with adjustments. This 
approach would allow an evaluation of the overall cost in- 
volved in solving the current two-dimensional problem. 


I - Win ” V’itf 
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The first solver is based on standard Gaussian elimination 
with partial pivoting and full storage. The second is based on 
triangular decomposition 8 and uses a compact storage scheme 
that avoids handling the zero entries and therefore should be 
more efficient than standard Gaussian elimination. The third 
solver is based on a Gauss-Scide! iterative scheme 9 and was not 
optimized for speed (through the choice of optimum accelera- 
tion parameters) but uses sparse matrix technology in proces- 
sing only the nonzero elements. The fourth and last solver used 
is based on the conjugate gradient method. 10 Handling the 
sparsity pattern for the third and fourth solvers is achieved by 
assembling the symbolic part of the coefficient matrix only 
once for a given grid size and given freestream (subsonic vs 
supersonic). The resultant structure is then stored on a disk 
file. Before the numerical part is executed, the symbolic infor- 
mation is read into the code and used directly to assemble the 
new matrix. This procedure is followed to reduce the time 
consumed in assembling the coefficient matrix. Notice also 
that in the Gauss-Seidel and conjugate-gradient solvers that 
the error tolerances for the coefficients involving maximum 
thickness, freestream Mach number, and location of maxi- 
mum camber were 1 .E-06, while those on angle of attack and 
maximum camber were l.E-04. 

Once the sensitivities of the potentials, and thus the Cp 
distribution, to the design variables are known, the sensitivity 
of the lift coefficients to the design variables can be easily 
computed. To minimize errors, these coefficients are com- 
puted using 

C t = 2T — 2((p u jz — v^ite) (20) 

and hence, 

dC L /dXD, = liS^jt/dXD, - d<p\jz/dXDj) (21) 

Finally, all methods used for computing the derivatives are 
compared to the finite-difference approach, and the results are 


Table 1 Accuracy of quasianalytical method 
for computing lift coefficient sensitivity derivatives 
for NACA 1406, GRID 81x20 


XD, 

Method 3 

? 

11 

O 

M 

M*. = 0.8 

Moo = 1 .2 

T 

FD 

0.0044 

0.5232 

-0.2949 


QA 

0.0044 

0.5447 

-0.3376 

M« 

FD 

0.0471 

0.9708 

1 .0235 


QA 

0.0470 

0.9905 

-0.0703 

a 

FD 

6.1385 

10.5229 

4.8758 


QA 

6.1386 

10.5229 

4.8726 

C 

FD 

9.9380 

19.5767 

0.7695 


QA 

9.9381 

18.6154 

0.7356 

L 

FD 

0.0696 

0.1499 

-0.0348 


QA 

0.0693 

0.1496 

-0.0349 


‘FD. finhe difference. QA. quasianalytical. 



x 

Fig. I Pressure coefficient distribution a! 0.2. 


presented and evaluated to determine the computational accu- 
racy and efficiency of each method. 

Test Cases 

In this sludy, the quasianalytical method has been used to 
determine the aerodynamic sensitivity coefficients at two 
freestream Mach numbers (M*, = 0.2, 0.8) for the NACA 1406 
airfoil at 1-deg angle of attack. Results were also obtained 411 
for a supersonic case at Af„= 1.2. Notice that further studies 
are needed to examine the results for a wider range of design 
parameter variation. 

In the following, two types of results will be presented. 
The first will be plots of Cp vs chord for the three cho- 
sen Mach numbers. The second will be the corresponding 

plots of (acp/ar>, (acp/aAC), (acp/a«), (acp/ac), and 

(dCp/dL) for the upper and lower surfaces and plots of 
(dACp/dT ), . . . , etc., involving the difference, all will be ob- 
tained by the quasianalytical method. In addition, all of the 
figures will also contain results obtained using the direct (finite 
difference) approach in which each design variable was indi- 
vidually perturbed by a small amount, typically 0.001, and a 
new flowfield solution obtained. Then the sensitivities were 
computed using ACp/AXD and are shown via dashed lines. In 
many cases, the lines are coincident with the quasianalytical 
results and cannot be observed. Table 1 compares results ob- 
tained by the two methods, and in most cases the agreement is 
within IVo. 

In all cases, an 81 x 20 stretched Cartesian grid was utilized. 
While finer grid studies are needed, they were not performed 
as part of this initial study. In addition, for these studies, the 
flowfield was normally computed using double precision arith- 
metic and the maximum residual reduced eight orders of mag- 
nitude. It was felt that this level of convergence was necessary 
in order to accurately evaluate sensitivity coefficients using a 
finite-difference approach, although such convergence may 
not be required in the flowfield solver for the quasi-analytical 
method. 
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Fig. 2 Sensitivity of pressure to maximum thickness, M r , 0.2. 
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Fig. 3 Sensilivity of pressure (o freestream Mach number. = 0.2. Flg ' 5 SensUivi,y of P ressurc <° maximum camber, Af» = 0.2. 




Results and Discussion 

Subsonic Case - =0.2 

Initial studies concentrated on subsonic cases since at least 
approximate results would be known from thin airfoil theory. 5 
Figure l shows the pressure distribution for the NACA 1406 
airfoil , while Figs. 2a and 2b show the sensitivity of the pres- 
sure to thickness for the same airfoil. As expected from thin 


airfoil theory, the upper and lower surface values are essen- 
tially identical, and the difference is very small everywhere. 
Also shown on the same figure (and on subsequent figures) by 
the dashed line is the result obtained by using the finite-differ- 
ence approach; and as can be seen, the agreement between the 
two approaches is excellent. 

The sensitivity of pressure to freestream Mach number is 
plotted on Figs. 3a and 3b. It is noticed that while the profiles 
for the upper and lower surfaces are similar, they are not equal 
in magnitude, indicating a nonlinear variation with Mach 
number as predicted by simple Prandtl-Glauert theory. How- 
ever, as indicated by the results plotted on Fig. 3b, the magni- 
tudes for this subsonic Mach number are very low. 

The sensitivity of the pressure coefficients to angle of attack 
are depicted for this case in Figs. 4a and 4b. As expected from 
linear thin airfoil theory, the upper and lower surface curves 
are essentially equal in magnitude but of opposite sign. Not 
surprisingly, the sensitivity of the delta Cp variation. Fig. 4b, 
has the shape of the pressure difference curve for a flat plate 
at angle of attack; and its magnitude, particularly near the 
leading edge, is quite large. 

On Figs. 5a and 5b is plotted the sensitivity of the pressure 
coefficient to the amount of maximum camber. Since camber 
contributes to lift, it is expected from the thin airfoil theory 
that these values should be “equal but opposite in sign’ 1 for 
the upper and lower surfaces. In addition, the pressure differ- 
ence curve has the correct shape for that associated with a 14 
mean line with the peak occurring at 30°/o chord 5 and has a 
magnitude comparable to those for the (dCp/da) curves. 

Finally, the sensitivity of pressure to the location of the 
maximum camber point is portrayed in Figs. 6a and 6b and, to 
say the least, the results are interesting. Since maximum cam- 
ber location affects the camber profile and hence lift, the equal 
and opposite behavior of the upper and lower surface coeffi- 
cients is expected. In addition, the pressure difference sensitiv- 
ity is primarily negative forward of the point of maximum 
camber and positive aft of it. This result indicates that if the 
location of maximum camber were moved rearward slightly 
(i.e., a positive A/.), that lift would be decreased on the for- 
ward portion of the airfoil and increased on the aft portion of 
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Fig. 6 Sensitivity of pressure to location of maximum camber. Fig. S Sensitivity of pressure to maximum thickness, Af<£ = 0.8. 

= 0 . 2 . 



Fig. 7 Pressure coefficient distribution at = 


the airfoil, which is in agreement with the results presented in 
Ref. 5. 

Transonic Case — = 0.8 

At A/« = 0.8 t the flow about the NACA 1406 airfoil has a 
strong shock at 40% chord, see Fig. 7, and the lower surface 
is entirely subcritical. As a consequence, the variation with 
chord of the sensitivity coefficients is considerably different 
than in the subsonic case. 

Figs. 8a and 8b show the sensitivity of pressure to the max- 
imum thickness; and while the lower surface profile is simi- 
lar to that obtained at subsonic conditions, the upper surface 
curve and the pressure difference coefficient plot show the 
effect of the upper surface shock wave. The large peak on the 
curves corresponds to the location of the shock wave and 
indicates that the shock-wave location is very sensitive to max- 
imum thickness. Notice on Figs. 8a and 8b the excellent agree- 
ment of the quasianalytieai results indicated by the solid lines 
with those obtained using the nnitc-diffcrcncc approach 
(dashed lines). 


The results for which are shown on Figs. 9a 

and 9b, are similar. The lower surface curve is typical of a 
subsonic flow, whereas the upper surface and the pressure 
difference coefficients reflect the presence of the upper surface 
shock wave. Similar comments can be made for the remaining 
design variable coefficients, which are plotted on Figs. 10, 11, 
and 12. 

Examination of the curves in the vicinity of the shock wave 
location indicates that the pressure sensitivity and indirectly 
the shock wave location is about equally influenced by the 
maximum thickness, freestream Mach number, and angle of 
attack. However, in comparison it is relatively insensitive to 
the location of maximum camber; but, perhaps surprisingly 
so, the pressure is twice as sensitive to the amount of maximum 
camber as it is to the other design variables. It should also be 
noticed that the lift is most sensitive to angle of attack and to 
maximum camber. 

In addition. Fig. 1 1 shows a discrepancy between the results 
obtained by the direct approach and those obtained through 
the quasianalytieai method. It will be shown in the follow- 
ing section that this discrepa ncy is related to the choice of the 
step size used in computing the finite-difference solution, thus, 
revealing a significant deficiency in computing the sensitiv- 
ity derivatives in nonlinear regimes via the finite-difference 
approach. 

Time Comparisons 

Obviously, in the development of the quasianalytieai 
method, it was hoped that not only would this approach yield 
accurate values for the aerodynamic sensitivity coefficients, 
but also that it would be more efficient than t he brute-force, 
finite-difference approach. Tabic 2 presents some comparisons 
concerning the amount of computational effort required to 
obtain solutions by the two approaches including results for 
the supersonic case. 4 M 

In comparing the values, several items should be kept in 
mind. First, it has been assumed that the finite-difference ap- 
proach will require six independent solutions. In practice, it 
might be possible to start each finite difference solution from 
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a previous solution and, thus, decrease the time to conver- 
gence. However, to be accurate, the finite-difference approach 
ms will probably require double precision and will have to be 
extremely well converged (i.e., l.E-08). Nevertheless, the val- 
ues for the finite-difference approach probably should be 
j|j viewed as maximum values. 

_ Second, the methods used for obtaining the sensitivity coef- 

ficients have not been optimized and, as mentioned earlier, 
may not even be optimum; and the flowfleld solution required 
tl for the quasianalytical approach may not need double preci- 
y sion and may not have to be as tightly converged. Thus, the 
values shown for the quasianalytical approach should also be 
viewed as maximum values. 

Us In spite of these limitations, results obtained by direct meth- 

y ods do indicate that the quasianalytical method is more com- 
putationally efficient at supersonic conditions and potentially 
efficient at transonic conditions than the brute-force, finite- 
difference approach. 

H| Notice that in this study, the initial guess used in computing 

the sensitivity derivatives via iterative methods was arbitrarily 
chosen as the zero vector. In addition, time comparisons pre- 
sented in Table 2 show that iterative methods arc in general less 


efficient than direct methods it the derivatives lor the current 
two-dimensional problem were sought about some general de- 
sign point. However, if the objective is to incorporate the 
sensitivity derivatives into an optimization loop (i.e., to use the 
derivatives in a continuation problem), then, a good initial 
guess (which in that case would be available) would enhance 
convergence, and the overall cost of computing the derivatives 
using iterative methods might be reduced. These points should 
be taken into consideration when a sensitivity study is to be 
integrated into an optimization procedure. 

Additional Test Cases 

The first group of cases are carried out to investigate the 
performance of the NACA 1406 airfoil for a range of Mach 
numbers from 0.79 to 0.86 in increments of 0.01 . As shown in 
Fig. 13. this range of transonic Mach numbers encompasses 
t lie development of the shock wave on the upper surface of the 
airfoil. Also, as shown on Figs. 14 and 15 for the cases involv- 
ing thickness, Mach number, and maximum camber, the 
quasianalytical derivatives arc in the vicinity of the shock wave 
frequently different from those obtained by the finilc-differ- 
cnce approach. This discrepancy raises two questions — what is 
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Tablt I Time* comparisons Tor obtaining 
sensitivity coefficients for five design variables 
for NACA 1406, GRID 81x20 


Method 6 

Moo = 0.2 

M* - 0.8 

M„ = 1.2 

FD 

1 .0000 

1.0000 

1.0000 

TD 

2.5187 

0.9962 

0.3929 

GE 

2. 4089 

0.9927 

0.5165 

GS 

0.9971 

1.5410 



CG 

35.2264 

10.6199 

— 


■All CPU limes were normalized by the lime taken to com- 
puic FD derivatives. 

b FD. finite difference; TD. triangular decomposition; 

GE, Gauss elimination; GS, Gauss-Seidcl; CG, conjugate 
gradient. 


Table 3 Effect of changing step sire delta on 
finite-difference lift coefficient sensitivity derivatives 
for NACA 1406, GRID 81x20, M<* = 0.84 


Delta XD, 

XD, = T 

XD, = Moo 

XD,=C 

I.E-03 

7.7603 

7.8715 

24.0912 

l.E-04 

-0.8493 

-0.6340 

83.9853 

I.E-05 

-0.8497 

-0.6364 

14.7719 

I.E-06 

-0.8498 

-0.6366 

14.7695 

QA a 

-0.8498 

-0.6367 

14.7692 


'QA, quaslanalytica! lift coefficient sensitivity derivatives. 
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Fig. 12 Sensitivity of pressure to location of maximum camber, 
- 0 . 8 . 
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Fig. 13 Pressure coefficient distributions at various Mach numbers. 
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Fig. 14 Finite-difference, pressure-difference coefficient sensitivity 
derivatives at various Mach numbers. 


the cause of the disagreement and which set of derivatives is 
more accurate? Examination of the variation of the integrated 
coefficient, dCL / d XD, with A/*, which is portrayed on Fig. 
16, shows that the quasianalyucal results arc smooth arid Tol- 
low a definite trend, whereas the finite-difference values are at 
best “discontinuous.” Consequently, it is concluded that the 
finite-difference results arc less accurate. 

In order to observe the performance of the finite-difference 
approach in the transonic regime, it is necessary to examine the 
effect of changing the stop size (delta of the design variable) on 
the computed derivatives, four different values for the step 
size (I.E-03, l.E-04, I.E-05, and l.E-06) were chosen and 


applied to the NACA 1406 at a Mach humBer of 0.84. Exam 
ination of this second group of results (see Tabic 3) show that 
as the step size is decreased, the finite-difference lift coefficient 
sensitivity derivatives approach the values computed by the 
quasianalytical method. However, in some cases, for small 
AXD t values, oscillations in the pressure coefficient sensitivity 
derivatives have been observed depending upon the machine 
used and the method of storing and retrieving the data. These 
oscillations combined with the difficulty of properly choosing 
a suitable finite-difference A XD, a priori indicates that the 
finite-difference approach is probably not a practical method 
of efficiently computing sensitivity coefficients. On the other 
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Fig. 15 Quasianalytical pressure-difference coefficient sensitivity de- 
rivatives at various Mach numbers. 



hand, the present results demonstrate that the quasianalytical 
method can be used accurately to obtain such coefficients in 
the transonic flight regime. 

Conclusion 

Based upon these investigations and results, it is concluded 
that the quasianalytical method is a feasible approach for ac- 
curately obtaining transonic aerodynamic sensitivity coeffi- 
cients in two dimensions. The results obtained from the quasi- 
analytical method are almost identical to those obtained by the 
brute-force (finite-difference) technique. Furthermore, the 
study indicates that the computation of sensitivity derivatives 
at transonic conditions is generally more accurate using the 
quasianalytical direct solver approach than the finite-differ- 
ence approach. In addition, the quasianalytical method is 
more efficient at supersonic Mach numbers and is potentially 
more efficient than the brute-force approach at transonic 
speeds. However, further studies to determine the effects of 
grid refinement and to examine the results over a wider range 
of design parameter variation are needed. 
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Abstract 

The quasianalytical approach is applied to the 
three-dimensional full potential equation to com- 
pute wing aerodynamic sensitivity coefficients in 
the transonic regime. Symbolic manipulation is 
used to reduce the effort associated with obtain- 
ing the sensitivity equations, and the large sen- 
sitivity system is solved using ” state of the art 
routines. Results are compared to those obtained 
by the direct finite difference approach and both 
methods are evaluated to determine their compu- 
tational accuracy and efficiency. The quasianalyt- 
ical approach is shown to be accurate and efficient 
for large aerodynamic systems. 

Nomenclature 

C Maximum camber in fraction of chord 

CG Conjugate gradient 

Cl Local lift coefficient 

CL Total lift coefficient 

Cp Pressure coefficient 

c(y) Chord function 

FD Finite difference 

GMRESGeneralized minimum residual 
L Chordwise location of maximum camber 

M Local Mach number 

M c Cutoff Mach number 0.94 < M c < l-° 

Moo Freestream Mach number 

Poo Freestream pressure, nondimensionalized 

by ( 27/(7 + 1)1 ft 
P 0 Stagnation pressure 

QA Quasianalytical 

goo Freestream velocity, nondimensionalized 

by V m 

T Maximum thickness in fraction of chord 

T 1..4 Twist angles 

U,V,W Contravariant velocity components 
V m Critical speed 

x,y,z Physical grid system 

X,Y,Z Computational coordinates 

* Graduate Research Asst. 

** Professor, Aerospace Engr., Associate Fellow AIAA 

Copyright c 1992 by the American Institute of 
Aeronautics and Astronautics, Inc. All rights re- 
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Leading edge function 

X-Coordinate of leading edge corner point 
X-Coordinate of trailing edge corner point 
Y-Coordinate of wing tip 
Vector of design variables 
Density, nondimensionalized by p 0 
Freestream density, nondimensionalized 
by po 

Stagnation density 
Retarded density coefficient 
First order backward difference operator 
Switching function 
Angle of attack 
Ratio of specific heats 
Reduced potential function 
$ Full potential function 

r Circulation 

Introduction 

To design transonic vehicles using optimiza- 
tion techniques requires aerodynamic sensitivity 
coefficients, which are defined as the derivatives 
of the aerodynamic functions with respect to the 
design variables. In most cases, the main con- 
tributor to the optimization effort is the calcu- 
lation of these derivatives; and, thus, it is de- 
sirable to have numerical methods which easily, 
efficiently, and accurately determine these coeffi- 
cients for large complex problems. At present 1 6 , 
there are two primary approaches for calculating 
transonic aerodynamic sensitivity derivatives. In 
the first approach, the sensitivities are calculated 
by perturbing a design variable from its previous 
value, a new complete solution is obtained, and the 
differences between the new and the old solutions 
are used to obtain the sensitivity derivatives. This 
brute force direct technique is computer intensive 
for complex governing equations that include a 
large number of design variables. In the second 
approach, termed the quasianalytical method, the 
sensitivities are obtained by solving a large sparse 
system of algebraic sensitivity equations in which 
the Jacobian matrix and right-hand-side vectors 
are obtained by differentiating the discretized form 
of the governing equations. The differentiations, 
while being staightforward in principle, are usually 
lengthy and tedious. However, once obtained, the 
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sensitivity equations can be very efficient and ac- 
curate for computing large numbers of sensitivity 
coefficients. 

In the first phase of this research 3 , the quasi- 
analytical approach was developed and applied 
to two-dimensional airfoils. Based upon these 
proof-of-concept investigations, it was concluded 
that the quasianalytical method was a feasible ap- 
proach for accurately obtaining transonic aerody- 
namic sensitivity derivatives in two dimensions, 
and was often more accurate and efficient than 
the finite difference method as the number of de- 
sign variables was increased. Further, the alge- 
braic forms of the matrix elements in the two- 
dimensional sensitivity equations were determined 
by hand, which involved extensive effort associated 
with differentiating the discretized residual with 
respect to the various design variables and the de- 
pendent unknowns. Today, such operations could 
be carried out using Symbolic Manipulation Pro- 
grams (SMs) 7 , such as MACSYMA* 1 ®, but present 
SMs are incapable of automatically performing 
all the necessary simplification, combinations, and 
cancellations of terms associated with algorithmic 
simplification of expressions. Consequently, the 
user must be familiar with the commands avail- 
able for the organization of expressions and con- 
duct various trials and experiments to identify a 
symbolic procedure which is efficient. As a result 
of these two-dimensional studies, it was decided 
to extend the quasianalytical approach to three 
dimensions and to investigate the use of Symbolic 
Manipulation Programs (SMs) l0,xl for obtaining 
the matrix elements. 


Problem Statement 

Application of the quasianalytical method to 
the full potential equation yields the sensitivity 
equation 

0R i,i '.k 1 ( 

^ j yaxD ) 1 ' 

where the residual expression in the computational 
plane in terms of backward differences is 

pU ? pV 

> pw 

+ ( 2 ) 

The retarded density coefficients in Eq.(2) are 





(3) 


(«) 

Pt,> t Jk + l/2 = + 1) 

(5) 


where 


Pij„k = [i - 7 (U + v *y + 7 . 1 ( 6 ) 

l 7 + 1 Ji 

and 


For this extended effort, it was decided 
to use for the flow solver a modified version 
of the three-dimensional direct-inverse analysis- 
design transonic full potential fully conservative 
code, ZEBRA 12 ” 15 . The full potential equation 
was selected because it can be solved rapidly and 
is robust, and accurate for engineering purposes 15 . 
Further, it can be formulated using a stretched 
Cartesian grid system that can be rapidly gener- 
ated and which has simple metrics. Also, such a 
grid permits the variation of several design param- 
eters without changing the physical or computa- 
tional grids. For the present work, the analysis 
portions of ZEBRA have been rearranged and un- 
needed portions deleted. In addition, the capabil- 
ity of calculating the sensitivity derivatives via the 
finite difference approach has been added. 
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In Eq.(7), the Mach number is obtained from 


PiJ 
and thus 


=d + 1 7 (8) 
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where Pij t k is nondimensionalized by p 0 . Prom 
Eqs.(7) and (9), 

Mij,k < i 
Mi,j,k > l 
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( 10 ) 



The contravariant velocities are 


U = (Xl + Xl)* x + Xy* Y (") 


p = fl - ^±(V* X + V*Y + W*z)\ ’ ‘ (22) 

l 7 + 1 J 


K = X v <tjr + *r ( 12 ) 

W = * z ( 1 3 ) 

where the full potential is split into perturbation 
and freestream components as 


and where the freestream values <joo, Px, and P, oo 
in Eqs.(20) and (21) are 


<7oo * I 
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Note that the angle of attack enters the formula- 
tion thru the above equation. Also note that the 
physical grid system (x,y,z) is transformed into the 
wing aligned computational grid (X,Y,Z) by 


X(x,y) 


x - zlc(y) 

c(y) 


(is) 


Y(y) = y 


(1C) 


Z(z) = z (17) 

The boundary conditions are the surface 
boundary condition, 

WmU fk +v $? (18) 

the Kutta condition along the wing semispan, 


r = XT '£ < x < oo (10) 

and the far field boundary condition. Additional 
conditions include updating the potential on the 
downstream boundary {<f>x = o) an( ^ implementing 
the wing symmetry condition by setting V = o. 

Once the unknown sensitivities d<j>/BXD are 
obtained, the sensitivities of the pressure coeffi- 
cient, Cp, with respect to the design variables can 
be computed. From the pressure coefficient ex- 
pression 


P - Poo 

Cp = — rrr 
wlor 1 


( 20 ) 


substitution for the pressure using the isentropic 
relation yields 


Cp = (7 + ^ (P 7 - plo) 

pqoo 


( 21 ) 


where 


Design variables can be classified according 
to whether or not they are coupled. Uncoupled 
design variables are termed basic variables, which 
are the independent variables that influence the 
solution of a problem; while coupled design vari- 
ables are termed nonbasic and are obtained from 
the basic design variables usually using simple al- 
gebraic expressions. For example, in the current 
problem, wing planform sweepback angles are non- 
basic design variables which are obtained knowing 
the basic variables or the coordinates of the cor- 
ner points of the wing. Other examples of nonbasic 
variables are the wing semi-span, aspect ratio, and 


taper ratio. 

The basic design variables for the current 


problem are: 

(a) Freestream design variables: These include 
the freestream Mach number and the angle of 
attack. The Mach number enters the formu- 
lation thru Eq.(23) while the angle of attack 
shows up in Eq.(14). 

(b) Cross section design variables: These include 
variables that define the airfoil section such as 
maximum thickness, maximum camber, and 
location of maximum camber for a NACA 
four-digit section and, variables that define 
each spanwise section such as geometric twist. 
For the current problem, these variables en- 
ter the problem via the boundary condition, 

Eq-(18). 

(c) Planform design variables: These variables 
define the geometry of the wing planform. In 
this study, the coordinates of the wing cor- 
ner points are used as the basic design vari- 
ables. Knowing the sensitivities with respect 
to these basic variables allows evaluation of 
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the derivatives with respect to the nonba- 
sic variables- The coordinates of the cor- 
ner points enter the current formulation via 
Eq.(15). 

Thus, for the current three-dimensional problem, 
the vector of design variables consists of twelve 
variables and six derived variables and is given by, 

XD — | Af , cl f T, C, L % T\ , T? t X L? , XT? t 

1 (2G) 

These variables are used in obtaining the right 
hand side vectors in Eq.(l). 

Symbolic and Numerical Treatment 

The basic approach used to symbolically dif- 
ferentiate the residual expression was to treat the 
main expression in terms of smaller subexpres- 
sions, each of which was examined in terms of 
its constituents. This process was extended until 
the final subexpressions included the appropriate 
derivative argument, the reduced potential or the 
design variables, in a simple functional form. The 
best method to obtain these subexpressions was to 
consider the governing equation and the involved 
intermediate expressions in the original form given 
in Eqs.(2)-(14). This splitting or nesting of expres- 
sions with various intermediate dependencies de- 
clared in advance allowed each subexpression to be 
handled efficiently by the symbolic manipulator. 
This usage of the chain rule of differentiation to- 
gether with MACSYMA’s ability to keep track of 
various equations resulted in an efficient scheme of 
analytical differentiation. It is noted that an early 
attempt to obtain the derivatives from a residual 
expressed as an explicit function of the reduced 
potential thru appropriate substitutions, Eq.(14) 
into (11), (12) and (13) up to Eq.(2), proved to 
be a poor strategy since the rapid increase in ex- 
pression size eventually caused MACSYMA to en- 
counter limitations on memory and manipulative 
ability. The experience gained from this attempt, 
however, turned out to be useful in identifying the 
capabilities and limitations of various MACSYMA 
commands and assisted in the development of fur- 
ther symbolic aspects associated with the project. 

During this study, various MACSYMA codes 
have been developed to assist in the application 
of the quasianalytical method. The first code, 
termed RMD.MAC, finds all residual reduced po- 
tential dependencies. This code is needed prior to 
carrying out the analytical differentiation of the 
residual, Eq.(2), with respect to the reduced po- 
tential function. Notice that the latter function 


shows up in Eq.(14), where the details of the de- 
pendence of the residual expression on this func- 
tion are not obvious, since intermediate expres- 
sions Eqs.(3) to (13) are involved. As mentioned 
earlier, handling each intermediate subexpression 
separately simplifies the operations involved. The 
result of this code is a file which includes various 
intermediate dependencies obtained in the form 
of lists. The second code termed RMDER.MAC, 
uses these lists and starts the symbolic differen- 
tiation process in order to obtain the Jacobian 
and right hand side vectors. The result of this 
lengthy code is a large FORTRAN segment that 
includes three subroutines and is about 15000 
lines long. As mentioned in the following section, 
this segment which is the heart of the quasian- 
alytical method, is linked into the quasianalyti- 
cal sensitivity driver. The third MACSYMA code 
is termed RCP.MAC, and generates FORTRAN 
source code for the derivatives of the pressure co- 
efficient, Eqs.(21) to (25), with respect to the vec- 
tor of design variables. This code uses the reduced 
potential sensitivity derivatives as input arrays. 
This segment of FORTRAN source code is also 
linked with the segment obtained from the second 
MACSYMA code. Finally, the fourth MACSYMA 
code is termed RESID.MAC and was created dur- 
ing debugging operations to test the evaluation of 
various residual terms. This program was very 
helpful in revealing logic and procedure errors in 
RMDER.MAC. Finally, it is important to empha- 
size that each of the above MACSYMA codes is 
executed only once followed by a transfer of the re- 
sulting FORTRAN segments to the QA sensitivity 
driver. 

Direct solvers that were previously used in the 
two-dimensional problem 2 (i.e. tridiagonal decom- 
position and full Gaussian elimination) failed on 
the three-dimensional problem due to limitations 
on memory; while the iterative routines developed 
earlier worked properly but were very slow. How- 
ever, library routines 16 available on the IBM-3090 
were extremely efficient with respect to memory 
and execution speed; and two scientific library 
solvers based on the iterative conjugate gradient 
method and the generalized minimum residual ap- 
proach have been used with success. For these 
solvers, the exact amount of storage needed de- 
pends on the sparsity and band width of the Ja- 
cobian matrix which in turn depends on the size 
of the three-dimensional grid. The present grid of 
45*30*16 yields a large, sparse, banded, and un- 
symmetric Jacobian matrix of about 17500*17500 
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that is less than one percent dense. An incomplete 
LU factorization is applied only once to this large 
matrix, and the sensitivity equations are solved 
using the iterative CG or GMRES methods 1 * 117 ’ 1 *. 
Following the factorization of the Jacobian matrix, 
back substitution using the known right hand side 
vectors generates the unknown sensitivity deriva- 
tives with a trivial computational cost. Recall that 
one crucial objective of this study is to exploit the 
efficiency of the QA method as the number of de- 
sign variables is increased. 

Program Structure 

The analysis-sensitivity program consists of 
the modified analysis program, ZEBRA, the finite 
difference sensitivity driver, and the quasianalyti- 
cal sensitivity driver. Execution of the main code 
starts with an analysis run followed by sensitivity 
derivative calculations for each point in the flow- 
field. These calculations are carried out either us- 
ing the FD method or the QA approach. The FD 
portion of the code uses two consecutive ZEBRA 
runs to calculate a vector of sensitivity derivatives. 
This brute force technique, while straight-forward, 
has the disadvantage of being expensive to im- 
plement and exhibits problems when single pre- 
cision variables are used. The QA driver consists 
of two main parts. The first part assembles the 
Jacobian matrix and the right-hand-side vectors. 
This assembly is achieved using calls to the large 
code segment generated via MACSYMA. This sec- 
tion of subroutines, as explained earlier, contains 
source code for the elements of the Jacobian ma- 
trix and right-hand-side vectors. Following the 
numerical assembly step, the second part of the 
sensitivity driver solves the sensitivity equations 
using one of the available linear sparse solvers and 
yields the unknown sensitivity vectors. Finally, 
the resulting sensitivity derivatives djtjdXD are 
processed to obtain the pressure coefficient sensi- 
tivity derivatives, dCpf dXD , at twenty-five chord- 
wise locations at each of the twenty wing semis- 
pan stations. This process is performed using the 
subroutines generated via RCP.MAC, the MAC- 
SYMA file used to symbolically differentiate the 
the pressure coefficient with respect to the reduced 
potential. 

Test Cases 

The wing configuration considered is that for 
the four cornered ONERA M6 wing planform 13 15 
with NACA 1406 airfoil sections. For this configu- 
ration, four test cases have been successfully con- 


ducted. The first case is subcritical at a freestream 
Mach number of 0.8 and an angle of attack of one 
degree. The second and third cases are supercrit- 
ical at Mach number of 0.84 and 0.88 respectively 
and an angle of attack of three degrees, while the 
fourth case is supersonic at a Mach number of 1.2 
and an angle of attack of three degrees. Due to 
space limitations, only results for the second case 
(Moo = 0 . 84 , a = 3 deg) are presented in this pa- 
per. This case is challenging since it includes a 
subcritical lower surface flow and exhibits an up- 
per surface shock wave located at 70% chord at 
the root to 10% chord at the tip that increases in 
strength from the root to a point near the wing 
tip. Thus, results for this case are believed suffi- 
cient to demonstrate the capabilities of the present 
analysis-sensitivity program. — 

In the above cases, a coarse-medium grid se- 
quence was used in computing the analysis infor- 
mation in order to speed up convergence. For 
the FD method, each design variable was indi- 
vidually perturbed by a small amount, typically 
1 * 10“*, and a new flowfield solution obtained. In 
all cases, double precision arithmetic was utilized 
and the residual reduced eight orders of magni- 
tude. In addition, the sensitivity information was 
computed by restarting each of the perturbed de- 
sign states from the coarse grid then proceeding 
to the medium grid. Different strategies for grid 
sequencing together with various choices of a suit- 
able starting solution are all valid options to speed 
up the FD approach. In the QA method, as men- 
tioned earlier, the sensitivity equation was set up 
with multiple right hand sides (the current vector 
of design variables, Eq.(26), includes twelve basic 
parameters) and was solved using the CG routine. 

Results and Discussion 

For the subcritical test case, the results ob- 
tained by the quasianalytical method were found 
to be in excellent agreement with results obtained 
via the finite difference method. In addition, the 
results followed the trend of the two-dimensional 
study 2 . 

Typical results for the chordwise variation of 
the pressure coefficient sensitivity derivatives at 
the Moo = 0 . 84 , a = 3 deg supercritjcal case are 
shown on Figs. I and 2 for a midspan station. 
Also displayed next to the legend in each case 
are the integrated coefficients dCi jbX As ex- 
pected, the sensitivity derivative profiles for the 
lower surface are typical of subcritical flow 2 ; and 
the upper surface results exhibit large variations 
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in the vicinity of the shock wave. The latter re- 
flect the influence on the aerodynamic coefficients 
of the sensitivity of the upper surface shock wave 
location to various design parameters. In addi- 
tion, a comparison of the spanwise distribution of 
the integrated coefficients, dCi/dX D tl is shown on 
Fig. 3; and in general these section values are 
smaller in magnitude than corresponding values 
for the two-dimensional problem 2 . As can be seen 
on Figs. 1-2, the agreement between the FD and 
QA results is excellent, indicating that accurate 
results can be obtained using the quasianalytical 
approch for three-dimensional transonic cases. 

Since both the FD and QA methods yield sim- 
ilar results, the question arises as to which is the 
least expensive. Current results obtained with the 
CG solver indicate that the QA method is com- 
putationally more efficient than the brute-force, 
finite difference approach. A representative CPU 
time ratio QA/FD of 0.46 was obtained for the 
above case with twelve design variables. It is rec- 
ognized, however, that the potential exists for re- 
ducing the cost associated with the FD method. 
For example, the perturbed runs could be exe- 
cuted directly on the medium grid starting with 
the design point solution obtained on the coarse 
grid. Likewise, the QA method could be improved 
by speeding up the evaluation of the Jacobian and 
right hand sides and/or using various options re- 
lated to the library solver. Therefore, the stated 
time ratio should only be considered as an esti- 
mate for comparing the two methods. 

One application of sensitivity derivatives is 
solution prediction. Fig. 4 compares the pres- 
sure coefficent distributions predicted at midspan 
using a first order Taylor series expansion about 
a design point with the actual pressure coefficient 
variation. The predicted Cp’s are calculated using 


dCp 


Cpprcdicted — C Pdctign + jj ^ XDi 


(27) 


where the dCp/dX D{ values where obtained from 
the QA approach and A XD{ = (0.005,0.1,0.005, 

0.001,0.1,0.1) for (Afoo , a, T t C,L,Ti) respectively. 

As can be seen, the agreement between the 
two distributions is very good, which indicates 
that the sensitivity derivatives calculated using 
the QA method can indeed be used in prediction. 
As mentioned earlier, another important applica- 
tion of sensitivity derivatives is in optimization 
routines. This application, however, is beyond the 
scope of the current project. 


Conclusion and Recommendations 

Based upon the above results, it is concluded 
that the quasianalytical method is a feasible and 
efficient approach for accurately obtaining tran- 
sonic aerodynamic sensitivity coefficients in three 
dimensions. In addition, use of the symbolic ma- 
nipulation package, MACSYMA, to carry out the 
symbolic evaluation of the elements of the sensi- 
tivity equations is crucial in this type of sensitiv- 
ity study. The results obtained from the quasi- 
analytical method are almost identical to those 
obtained by the finite difference approach. Fur- 
thermore, the study indicates that (a) obtaining 
the quasianalytical sensitivity derivatives using an 
iterative conjugate gradient method is more ef- 
ficient than computing the derivatives by the fi- 
nite difference method, especially as the number 
of design variables increases, and (b) the quasi- 
analytical method shows promise with regard to 
analysis-sensitivity methodologies applied to large 
aerodynamic systems. 

Acknowledgment 

This project is supported by the National 
Aeronautics and Space Administration under 
grant No. NAG 1-793, Dr. E. Carson Yates, Jr., 
and Dr. Woodrow Whitlow, Jr., NASA Langley 
Research Center, as technical monitors. The au- 
thors also appreciate the helpful suggestions and 
discussions of Dr. Perry Newman of NASA Lang- 
ley. 

References 

1. Sobieski, J.S. "The Case for Aerodynamic 
Sensitivity Analysis”, Paper presented to 
NASA/ VPI&SU Symposium on Sensitivity 
Analysis in Engineering, September 25-26, 
1986. 

2. Elbanna, H.M. and Carlson, L.A., "Deter- 
mination of Aerodynamic Sensitivity Coeffi- 
cients Based on the Transonic Small Pertur- 
bation Formulation”, Journal of Aircraft , Vol 
27, No. 6, June 1990, pp 507-515. 

3. Bay sal, 0., Eleshaky, M.E., Burgreen, G.W., 
"Aerodynamic Design Optimization Using 
Sensitivity Analysis and Computational Fluid 
Dynamics”, AIAA paper No. 91-0471. 

4. Dulikravich, G.S.,” Aerodynamic Shape De- 
sign and Optimization”, AIAA paper No. 91- 
0476. 

5. Taylor III, A.C., Korivi, V.M., Hou, G.W., 
"Approximate Analysis and Sensitivity Anal- 


6 



ysis Methods for Viscous Flow Involving Vari- 
ation of Geometric Shape”, AIAA paper 
No. 91-1569. 

6. Bay sal, 0., Eleshaky, M.E., Burgreen, G.W., 

” Aerodynamic Shape Optimization Using 
Sensitivity Analysis on Third-Order Euler 
Equations”, AIAA paper No. 91-1577. 

7. Bau, H.H., ” Symbolic Computation - An 
Introduction for the Uninitiated , Symbolic 
Computation in Fluid Mechanics and Heat 
Transfer, The American Society of Mechan- 
ical Engineers, Vol 97, 1988, pp 1-10. 

8. Macsyma Reference Manual, Version 13, 
Computer Aided Mathematics Group, Sym- 
bolics, Inc., 1988. 

9. Macsyma User’s Guide, Computer Aided 
Mathematics Group, Symbolics, Inc., 1988. 

10. Roach P. and Steinberg S., "Symbolic Manip- 
ulation and Computational Fluid Dynamics” , 
AIAA Journal , Vol 22, No. 10, 1984, pp 1390- 
1394. 

11. Steinberg S. and Roach P., "Automatic Gen- 
eration of Finite Difference Code”, Symbolic 
Computation in Fluid Mechanics and Heat 
Transfer, The American Society of Mechan- 
ical Engineers, Vol 97, 1988, pp 81-86. 

12. South Jr., J.C., Keller, J.D., Hafez, M.M., 
"Vector Processor Algorithms for Transonic 
Flow Calculations”, AIAA Journal , Vol 18, 


No. 7, 1980, pp 786-792. 

13. Weed, R.A., Anderson, W.K., Carlson, L.A., 
"A Direct-Inverse Three Dimensional Tran- 
sonic Wing Design Method for Vector Com- 
puters”, AIAA paper No. 84-2 156. 

14. Carlson, L.A. and Weed, R.A., ”A Di- 
rect Inverse Transonic Wing Design Analysis 
Method with Viscous Interaction” , AIAA pa- 
per No. 85-4075. 

15. Carlson, L.A. and Weed, R.A.,^ "Direct- 
Inverse Transonic Wing Analysis Design 
Method with Viscous Interaction”, Journal of 
Aircraft , Vol 23, No.9, Sept.1986, pp 711-718. 

16. IBM Engineering and Scientific Subroutine 
Library, Guide and Reference, Rel 3, SC23- 
0184-3. 

17. Saad, Y. and Schultz, M.H., ”GMRES: A 
Generalized Minimum Residual Algorithm 
for Solving Nonsymmetric Linear Systems”, 
SIAM Journal of Scientific and Statistical 
Computing , Vol 7, No.3, 1986, pp 856-869. 

18. Sonneweld, Wesseling, and De Zeeuv, Multi- 
grid and Conjugate Gradient Methods as 
Convergence Acceleration Techniques in Mul- 
tigrid Methods for Integral and Differential 
Equations, pp 117-167, Edited by Paddon, 
D.J. and Holstein, M., Oxford University 
Press (Claredon), Oxford. 


7 






















Fig. 2 Pressure Coefficient Sensitivity Derivatives at 56 % Span 




Fig. 3 Sensitivity of Cl to M and a 
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Fig. 4 Actual vs Predicted Pressure Coefficients at 56 % Span 
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M 

The quasianalytical approach is applied to the three-dimensional full potential equation 
l , to compute wing aerodynamic sensitivity coefficients in the transonic regime. Symbolic ma- 

I nipulation is used to reduce the effort associated with obtaining the sensitivity equations, 

and the large sensitivity system is solved using ”state of the art” routines. The quasian- 
alytical approach is believed to be reasonably accurate and computationally efficient for 
H three-dimensional problems. 


INTRODUCTION 


To design transonic vehicles using codes which utilize optimization techniques requires 
aerodynamic sensitivity coefficients, which are defined as the derivatives of the aerodynamic 
functions with respect to the design variables. In most cases, the main contributor to the 
optimization effort is the calculation of these derivatives; and, thus, it is desirable to have 
numerical methods which easily, efficiently, and accurately determine these coefficients for 
large complex problems. The primary purpose of the present study is to investigate the 
application of the quasianalytical method [1,2] to three-dimensional transonic flows using as 
the fundamental flow solver the three-dimensional transonic full potential fully conservative 
code, ZEBRA [3]. 


PROBLEM STATEMENT 


Application of the quasianalytical method to the full potential equation yields the sen- 
sitivity equation 


dRij t k ( _ ( dR iJik \ 

afe.fc* \ ax d ) “ ^ axo ) 


0 ) 



where XD is the vector of design variables and the residual expression, of the full 

potential equation in the computational plane, X,Y, Z, in terms of backward differences is 


c 
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Ri,j,k — &X (^-7—); + 1/2.J ,fc + &Y[—j-)i,}+ 1/2.1 + &z( — )i,j,k+ 1/2 


c ,PW, 
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( 2 ) 


Here, the retarded density p and the contravariant velocity components U,V, and W, are 
lengthy functions of the reduced potential function, <j>. The boundary conditions for Eq.(2) 
are the surface condition, W = w + the Kutta condition along the wing semispan, 

T = A <f>, x TE < x < 00, and the fartield condition. Additional conditions are the downstream 
boundary potential 4> x = 0 and the wing symmetry condition, V = 0. 

The discretized form of Eq.(2) contains lengthy expressions, and mathematical symbolic 
manipulation [4-6] was used to determine the functional dependencies of the residual, the 
analytical forms of the derivatives, and to generate the corresponding computer code. The 
basic approach used to differentiate the residual expression was to treat the main expression 
in terms of smaller subexpressions, each of which was examined in terms of its constituents. 
This process was extended until simple functional forms for the derivatives were obtained. 
This subdivision and chain rule differentiation by symbolic manipulation efficiently generated 
source code for the jacobian and vectors in Eq.(l). The resultant large sparse system, 
typically 17500 * 17500, of algebraic equations is then efficiently solved for 5^ using either 
the iterative conjugate gradient method or the generalized minimum residual algorithm [7- 
8]. From these, the pressure and lift coefficient sensitivities to the design variables can be 
computed. Notice that the effort associated with this approach is essentially independent of 
the number of design variables considered on the right-hand-side of Eq.(l). 


EXAMPLE AND DISCUSSION 


Consider the ONERA M6 wing planform with NACA 1406 airfoil sections at a super- 
critical condition of Af TO = 0.84 and a = 3 deg, which has sub critical lower surface flow and 
exhibits an upper surface shock wave located at 70 % chord at the root to 10 % chord at 
the tip that increases in strength from the root to a point near the wing tip. Basic design 
variables for the current problem include freestream design variables, Mach number M x and 
angle of attack a; cross-section design variables of maximum thickness, T, maximum camber, 
C, and location of maximum camber, L] variables that define wing twist, T1.T2.T3, and T4; 
and planform tip coordinates, X LE tip , XTE tip , and Y tip . Knowing the sensitivities to these 
basic design parameters permits subsequent evaluation of the derivatives with respect to 
the nonbasic variables taper ratio, aspect ratio, wing area, and sweepback angles. Thus, the 
present method determines sensitivity coefficients for twelve design variables and five derived 
design variables. 

As part of the solution d<f>/dXD values are obtained for every grid point in the flowfield. 
Also, the method automatically computes dC p /dXD at twenty-five chord wise locations at 
each of the twenty semi-span stations on the wing as well as dQ/dXD at each of the span 
stations. Typical results for the example case are shown in Fig.l for a midspan station. As 
expected, the sensitivity derivative profiles for the lower surface are typical of subcritical 
flow [2]; and the upper surface results exhibit large variations in the vicinity of the shock 
wave. The latter reflect the influence on the aerodynamic coefficients of the sensitivity of 
the upper surface shock wave location to the various design parameters. Currently, efforts 
arc in progress to validate the present method by comparison with the finite difference 
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approach, which calculates sensitivities by perturbing a design variable from its previous 
value, obtaining a new solution, using the differences between the new and old solutions 
to obtain the sensitivity coefficients. While this direct technique is computer intensive and 
inefficient, it should serve as a check on the present method. 

Based upon the present results, it is concluded that the quasianalytical method is a 
viable and efficient concept for the determination of three-dimensional transonic aerodynamic 
sensitivity coefficients. In addition, use of symbolic manipulation to evaluate the elements 
of the sensitivity equation is believed to be an efficient approach to the development of 
such methods. Finally, further studies are needed to determine the accuracy and range of 
applicability of the quasianalytical approach. 
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Abstract 

The quasi-analytica! approach is applied to the three-dimensional full potential equation to compute wing 
aerodynamic sensitivity coefficients in the transonic regime. Symbolic manipulation is used and is crucial in 
reducing the effort associated with obtaining sensitivity equations, and the large sensitivity system is solved using 
sparse solver routines such as the iterative conjugate gradient method. The results obtained are almost identical to 
those obtained by the finite difference approach and indicate that obtaining the sensitivity derivatives using the quasi- 
analytical approach Is more efficient than computing the derivatives by the finite difference method, especially as 
the number of design variables increases. It is concluded that the quasi -analytical method is an efficient and accurate 
approach for obtaining transonic aerodynamic sensitivity coefficients in three dimensions. 
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Introduction 

To design transonic vehicles using optimization techniques requires aerodynamic sensitivity coefficients, 
which are defined as the derivatives of the aerodynamic functions with respect to the design variables. In most cases, 
the main contributor to the optimization effort is the calculation of these derivatives; and, thus, it is desirable to have 
numerical methods which easily, efficiently, and accurately determine these coefficients for large complex problems. 
At present* there are two primary approaches for calculating transonic aerodynamic sensitivity derivatives. In the 
first approach, the sensitivities are calculated by perturbing a design variable from its previous value, a new 
complete solution is obtained, and the differences between the new and the old solutions are used to obtain the 
sensitivity derivatives. This brute force direct technique is computer intensive for complex governing equations that 
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include a large number of design variables. In the second approach, termed the quasi-analytical method, the 
sensitivities are obtained by solving a large sparse system of algebraic sensitivity equations. While the matrix 
elements in these algebraic sensitivity equations are obtained analytically, they are obtained by analytically 
differentiating the discretized or numerical forms of the equations governing the flowfield. Further, the aerodynamic 
and sensitivity solutions are obtained numerically. Thus, the method is termed a quasi-analytical rather than a 
numerical or analytical method. It should be noted that the differentiations to obtain the coefficients for the 
algebraic sensitivity equations, while being straightforward in principle, are usually lengthy and tedious. However, 
once obtained, the sensitivity equations can be very efficient and accurate for computing large numbers of sensitivity 
coefficients. 

In the first phase of this research 2 , the quasi-analytical approach was developed and applied to two- 
dimensional airfoils. Based upon these proof-of-concept investigations, it was concluded that the quasi-analytical 
method was a feasible approach for accurately obtaining transonic aerodynamic sensitivity derivatives in two 
dimensions and was often more accurate and efficient than the finite difference method as the number of design 
variables was increased. Further, the algebraic forms of the matrix elements in the two-dimensional sensitivity 
equations were determined by hand, which involved extensive effort associated with differentiating the discretized 
residual with respect to the various design variables and the dependent unknowns. Today, such operations could be 
carried out using symbolic manipulation programs^, such as MACSYMA 1 ^’ 11 , but present symbolic manipulators 
are incapable of automatically performing ail the necessary simplification, combinations, and cancellations of terms 
associated with algorithmic simplification of expressions. Consequently the user must be familiar with the commands 
available for the organization of expressions and conduct various trials and experiments to identify a symbolic 
procedure which is efficient. As a result of these two-dimensional studies, it was decided to continue the research. 
Consequently, the primary objectives of the present effort have been to apply the quasi-analytical method to three- 
dimensional transonic flow, investigate the use of symbolic manipulation programs 12,13 for obtaining the matrix 
elements of the sensitivity equations, and to determine the efficiency and accuracy of the quasi-analytical approach 
for determining transonic aerodynamic sensitivities. 

For this extended effort, it was decided to use for the flow solver a modified version of the three- 
dimensional direct-inverse analysis-design transonic full potential fully conservative code, ZEBRAII 14 ' 17 . The full 
potential equation was selected because it can be solved rapidly and is robust and accurate for engineering 
purposes 17 . Further, it can be formulated using a stretched Cartesian grid system that can be rapidly generated and 
which has simple metrics. Also, such a grid permits the variation of several design parameters without changing 
the physical or computational grids. For the present work, the analysis portions of ZEBRAII have been rearranged 
and unneeded portions deleted. In addition, the capability of calculating the sensitivity derivatives via the finite 
difference approach has been added. 
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Problem Statement 


Application of the quasi-analytical method to the full potential equation yields the sensitivity equation 

(i) 


dR . .. 
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where the residual expression in the computational plane in terms of backward differences is 
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The retarded density coefficients in Eq. (2) are 
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In Eq. (7), the Mach number is obtained from 
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and thus 
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where p is nondimensionalized by pQ. From Eqs. (7) and (9), 
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The contravariant velocities are 




( 11 ) 


+ < *r 


w= <x>. 


( 12 ) 

(13) 


and the full potential is split into perturbation and freestream components as 

+ Xq. Cos ( a ) + Zq_Sin(a) 


(14) 


Note that the angle of attack enters the formulation thru the above equation and that the physical grid system (x,y,z,) 
is transformed into the wing aligned computational grid (X,Y,Z) by 


X(x>y) = £^M 

c(y) 


(15) 


Y(y)=y 


(16) 
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Z(z) = z 


(17) 


The boundary conditions are the surface boundary condition, 


W = U— + V— 
dX dY 


(18) 


the Kutta condition along the wing semispan. 


T = A<J), x TE <x<. 


(19) 


and the farfield boundary condition. Additional conditions include updating the potential on the downstream 
boundary ( <f> x =0) and implementing the wing symmetry condition by setting V = 0. 

Once the unknown sensitivities 3<J >/3XD are obtained, the sensitivities of the pressure coefficient, Cp, with 
respect to the design variables can be computed. From the pressure coefficient expression 


Cp = 


P-P m 

P qlP- 


( 20 ) 


substitution for the pressure using the isentropic relations yields, 
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where p is given by Eq. (6) and where the freestream values and P m in Eqs. (20) and (21) are 
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Design Variables 

Design variables can be arbitrarily classified according to whether or not they are coupled. Uncoupled 
design variable are termed basic variables and are the independent variables that influence the solution of a problem. 
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Coupled design variables are defined here to be nonbasic and are obtained from the basic design variables usually 
using simple algebraic expressions. For example, in the current problem, wing planform sweepback angles are 
defined as nonbasic design variables since they are obtainable from the basic variables, i.e. the coordinates of the 
comer points of the wing. Other examples of nonbasic design variables are the wing area, aspect ratio, and taper 
ratio. 

The basic design variables for the current problem include freestream variables, airfoil cross-section 
variables, and planform parameters. The freestream design variables include the freestream Mach number, M aa , 
and the angle of attack, a. The Mach number enters the formulations thru Eq. (22) while the angle of attack shows 
up in Eq. (14). The airfoil section design variables include maximum thickness, maximum camber, location of 
maximum camber, and four angles that define at each spanwise station the amount of geometric twist. These 
variables enter the problem via the wing surface boundary condition, Eq. (18). The basic planform design variables 
define the geometry of the wing planform and are comprised of the coordinates of the wing corner points, which 
enter the formulation via Eq. (15). Evaluation of the sensitivities with respect to these basic planform variables 
allows the determination of the derivatives with respect to the nonbasic variables. Thus, for the current three- 
dimensional problem, the vector of design variables consists of twelve basic variables and is given by, 

XD = a, T, C,L,T V T z , T v T 4 ,XL t ,XT t , ( 25 ) 

These variables are used in obtaining the right hand side vectors in Eq. (1). 

Note that the design variables listed in Eq. (25) form a complete set of the basic design variables 
influencing the aerodynamic solution for the wing planform and wing sections considered in the present 
investigation. If the wing planform were more complicated, having for example a leading edge break, then the 
coordinates of that break point would have to be included in the vector of design variables. For more complex 
configurations or for problems involving coupling such as aeroelastic phenomena, the design variable set would be 
found by examining the solution model(s) and determining which flowfield and geometric parameters appear and 
consequently affect the aerodynamic solution. 

Symbolic and Numerical Treatment 

The basic approach used to symbolically differentiate the residual expression was to treat the main 
expression in terms of smaller subexpressions, each of which was examined in terms of its constituents. This process 
was extended until the final subexpressions included the appropriate derivative argument, the reduced potential or 
the design variables, in a simple functional form. The best method to obtain these subexpressions was to consider 
the governing equation and the involved intermediate expressions in the original form given in Eqs. (2)-(I4). This 
splitting or nesting of expressions with various intermediate dependencies declared in advance allowed each 
subexpression to be handled efficiently by the symbolic manipulator, in this case MACSYMA. This usage of the 
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chain rule of differentiation together with the ability of the symbolic manipulator to keep track of various equations 
resulted in an efficient scheme of analytical differentiation. It is noted that an earlier attempt to obtain the derivatives 
from a residual expressed as an explicit function of reduced potential thru appropriate substitutions, Eq. (14) into 
(11), (12) and (13) up to Eq. (2), proved to be a poor strategy since the rapid increase in expression size eventually 
caused the manipulator program to encounter limitations on memory and manipulative ability. The experience gained 
from this attempt, however, turned out to be useful in identifying the capabilities and limitations of various symbolic 
commands and assisted in the development of further symbolic aspects associated with the project. 

During this study, various symbolic manipulator codes were developed to assist in the application of the 
quasi-analytical method. The first code, found all residual reduced potential dependencies. This code was needed 
prior to carrying out the analytical differentiation of the residual, Eq.(2), with respect to the reduced potential 
function. Notice that the latter function shows up in Eq. (14), where the details of the dependence of the residual 
expression on this function are not obvious, since intermediate expressions Eqs.(3) to (13) are involved. As 
mentioned earlier, handling each intermediate subexpression separately simplifies the operations involved. The result 
of this code was a file which included various intermediate dependencies obtained in the form of lists. The second 
code used these lists to perform the symbolic differentiation process to obtain the Jacobian and right hand side 
vectors for Eq. (1), and the result of this lengthy code was a large 15000 line FORTRAN segment that included 
three subroutines. This segment is the heart of the quasi-analytical method and is linked into the quasi-analytical 
sensitivity driver. The third symbolic code generated FORTRAN source code for the derivatives of the pressure 
coefficient, Eqs. (21) to (25), with respect to the vector of design variables and used the reduced potential sensitivity 
derivatives as input arrays. This segment of FORTRAN source code was then also linked with the segment obtained 
from the second symbolic code. Finally, the fourth symbolic code was created during debugging operations to test 
the evaluation of various residual terms and was very helpful in revealing logic and procedure errors. Finally, it 
is important to emphasize that each of the above symbolic codes is executed only once followed by a transfer of 
the resulting source segments to the quasi-analytical sensitivity driver. Details and sample MACSYMA codes for 
these processes are given in Ref. 4. 

Direct solvers that were previously used in the two-dimensional problem 2 failed on the three-dimensional 
problem due to limitations on memory; while the iterative routines developed earlier worked properly but were very 
slow. However, library solvers 18 based on the iterative conjugate gradient method and the generalized minimum 
residual approach have been used with success and have proven to be extremely efficient with respect to memory 
and execution speed. For these solvers, the exact amount of storage needed depends on the sparsity and band width 
of the Jacobian matrix which in turn depends on the size of the three-dimensional grid. The present grid of 45 x 
30 x 16 yields a large, sparse, banded, and unsymmetric Jacobian matrix of (43 X 29 X 14) X (43 X 29 X 14) or 
about 17500 x 17500 that is less than one percent dense. An incomplete LU factorization is applied only once to 
this large matrix, and the sensitivity equations are solved using the iterative CG or GMRES methods 18 ’ 19 ’ 20 . 
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Following the factorization of the Jacobian matrix, back substitution using the known right hand side vectors 
generates the unknown sensitivity derivatives with a trivial computational cost. This approach exploits the efficiency 
of the QA method as the number of design variables is increased. 

Program Structure 

The analysis-sensitivity program consists of the modified flowfield analysis program, ZEBRA, the finite 
difference sensitivity driver, and the quasi-analytical sensitivity driver. Execution of the main code starts with an 
analysis run followed by sensitivity derivative calculations carried out either using the FD method or the QA 
approach. The FD portion of code uses two consecutive ZEBRA runs to calculate a vector of sensitivity derivatives. 
This brute force technique, while straight-forward, has the disadvantage of being expensive to implement and 
exhibits problems when single precision variables are used. The QA driver consists of two main parts. The first part 
assembles the Jacobian matrix and the right -hand -side vectors thru calls to the large code segment generated via 
symbolic manipulation. This section of subroutines, as explained earlier, contains source code for the elements of 
the Jacobian matrix and right-hand-side vectors. Following the numerical assembly step, the second part of the 
sensitivity driver solves the sensitivity equations using one of the available linear sparse solvers and yields the 
unknown sensitivity vectors at each point in the flowfield. Finally, the resulting sensitivity derivatives, dfyjdXD , 
are processed to obtain the pressure coefficient sensitivity derivatives SCpfdXD , at twenty-five chordwise locations 
at each of the twenty wing semispan stations. 

Test Cases 

For the present study, most of the test cases utilized the four cornered ONERA M6 wing planform 1517 
with a variety of airfoil sections including NACA 1406, 1706, 2406, and 2706 airfoils. This planform has an aspect 
ratio of 3.8 and a taper ratio of 0.56, with leading and trailing edge sweeps of thirty and 15.76 degrees respectively. 
Freestream conditions included subcritical cases at Mach 0.8 and an angle of attack of one degree, several 
supercritical transonic cases, and some supersonic cases up to Mach 1.2. Due to space limitations, most of the 
results of this paper will be for the ONERA M6 planform with NACA 1406 airfoil sections at freestream conditions 
of Mach 0.84 and three degrees angle of attack. This case is challenging since it has a subcritical lower surface flow 
and exhibits an upper surface shock wave located at 70% chord at the root that shifts to 10% chord at the tip and 
which increases in strength from the root to a point near the wing tip. Thus, the results for this case should be 
sufficient to demonstrate the capabilities of the present analysis-sensitivity method at transonic conditions. Complete 
detailed results for all the cases are presented in Refs. 4 and 21. 

In the above cases, a coarse-medium grid sequencing was used in the flowfield computations to enhance 
convergence. For the finite difference method of computing the sensitivities, each design variable was individually 
perturbed a small amount, typically 1 x lCT 6 , and a new flowfield solution obtained. In all cases, double precision 
arithmetic was utilized and the residual reduced eight orders of magnitude. In addition, the finite difference 
sensitivity results were computed by restarting each of the perturbed design states from the coarse grid then 
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proceeding to the medium grid. Different strategies for grid sequencing, such as starting on the medium grid with 
a previously obtained converged soulution, are all valid options to speed up the finite difference approach; but these 
were not investigated in this study. In the quasi-analytical method, the sensitivity equation, Eq. (1), was solved with 
twelve right hand sides representing the vector of design variables, Eq. (25), using one of the sparse solvers. In all 
cases, dfyjdXD values were obtained for every gridpoint in the flowfield. Also, the method automatically computed 
upper and lower surface dCp/dXD values at twenty-five chordwise locations at each of the twenty semispan stations 
on the wing planform as well as the dCljdXD values at each of the span stations and the overall win gdCL/dXD 
results. 

Results and Discussion 

For the subcritical test cases, the results obtained by the quasi-analytical method were found to be in 
excellent agreement with results obtained from the finite difference method. In addition, the results followed the 
trend of the two-dimensional study. 2 

Representative results for the chordwise variation of the pressure coefficient and its sensitivity derivatives 
for a supercritical case (M^ = 0.84, or = 3°) are shown on Fig. 1 for 56.4 percent semispan. Displayed in the 
comer in each case are the integrated coefficients, section lift coefficient for the pressure distribution and dCl/dXD i 

for the rest. In subcritical flow, the sensitivities with respect to the Mach number and the thickness would be small 
and similar for the upper and lower surfaces while those for angle of attack, camber, and camber location would 
have larger upper and lower surface values of opposite sign. As expected the sensitivity derivative profiles on Fig. 

1 for the lower surface are typical of subcritical flow 2 . However, the upper surface results exhibit large variations 
in the vicinity of the shock wave that reflect the influence on the pressure of the sensitivity of the upper surface 
shock wave location to various design parameters. As can be seen by noting the differences in vertical scale, the 
pressures and lift coefficient at this mid semispan location are most sensitive to camber and least sensitive to camber 
location. Finally, the agreement between the finite-difference and quasi-analytical predictions is excellent, indicating 
that accurate three-dimensional transonic results can be obtained using the quasi-analytical approach. 

Some results for the spanwise variation of the section lift sensitivity derivatives are shown on Fig. 2, where 
the numbers in the lower left comer in this case are the total wing lift coefficient sensitivities, dCLjdXD Note that 
the sensitivity of section lift to freestream Mach number and angle of attack is relatively constant over most of the 
semispan, but that the lift sensitivity to wing twist at twenty (T2) and sixty (T3) percent semi -span is concentrated 
in the region near the twist location. While not shown, lift sensitivities to twist at the root and the wing tip are only 
one-third to one-fourth of those at midspan. In general, primarily due to wing sweep and finite span, all the section 
sensitivity values are smaller in magnitude than corresponding values for the two-dimensional problem 2 . Finally the 
agreement between the finite difference and quasi-analytical section sensitivities is excellent. 

Figure 3 shows representative section and wing lift sensitivity derivative results for some of the nonbasic 
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design variables. While the total wing lift sensitivities can often be obtained by other means, the present method 
also yields spanwise and chordwise information. Note that while the lifts are relatively insensitive to the semispan, 
the outboard lift and total lift exhibit a strong dependence on area, aspect ratio, and taper ratio, and that the 
agreement between the quasi-analytical method and the finite difference approach is reasonable. While not shown, 
the corresponding derivatives with respect to the leading and trailing edge sweep angles were very small. 4 » 21 

While both the finite difference and quasi-analytical methods yield similar results, the present results 
indicate that for twelve design variables the quasi-analytical method is about 2.4 times computationally more efficient 
than the brute force finite difference approach. However, it is recognized that the costs associated with the finite 
difference method probably could be reduced by executing the perturbed runs directly on the medium grid starting 
with the design point solution obtained from the coarse grid. Likewise, the quasi-analytical method could be 
improved by utilizing various options associated with the sparse system equation solvers; and both methods are 
probably affected by grid size. Therefore, the stated relative efficiency should only be viewed as an estimate when 
comparing the two methods. 

One application of sensitivity derivatives is solution prediction, and Fig. 4 compares two pressure 
coefficient distributions at the 56 percent semispan location predicted using a first order Taylor series expansion 
about the original calculation point with the actual variation. The predicted Op’s were calculated using 


Cp pM = CPoriginai + (dCpfdXD) A XD t 


( 26 ) 


where the dCpjdXD i values were obtained from the quasi-analytical method at Mach 0.84 and a = 3 degrees. 

For the two results presented, the wing thickness was increased 8.3 percent to 0.065 chord and the wing tip leading 
edge ordinate was moved aft 0. 1 chord respectively. Since the original lift coefficient at this station, as shown on 
Figure 1(a), was 0.383, both changes resulted in a slight increase in lift coefficient and aft movement of the 
shockwave at this station. However, the detailed results 4,21 show that the movement of the wing tip ordinate caused 
a lift coefficient decrease in the inboard sections of the wing. As can be seen on the figure, the agreement between 
the quasi -analytically predicted and actual pressure distributions is very good, which indicates that the sensitivity 
derivatives calculated using the quasi-analytical method can be used for predictions. Similar results were obtained 
for the other design variables. 4 

Since sensitivity derivatives describe the response of the overall solution to changes in design variables, 
they can be computed over a range of flight conditions to determine the degree and nature of the influence of each 
design variable on the solution. At transonic conditions, the Mach number strongly influences a wing fiowfield; and, 
thus, sensitivity derivatives were computed for the ONERA M6 planform with NACA 1406 cross sections at an 
angle of attack of three degrees for Mach numbers ranging from 0.8 to 1.2. For simplicity, only the derivatives of 
the total wing lift coefficients with respect to each of the twelve basic design variables were considered. Figure 5 




shows results for three of these design variables, freestream Mach number, maximum camber location, and wing 
tip trailing edge ordinate. For all the design variables the largest variation of each derivative occurs in the transonic 
regime below Mach one. In this range as Mach number increases, the upper surface shock wave is rapidly moving 
towards the trailing edge with the inboard portions reaching the trailing edge, first. Thus, as shown by dCLjdM^ 
there initially is an increase in wing lift coefficient. However, by Mach 0.92, the inboard portion of the shock wave 
is at or near the trailing edge, and the effects of lower surface pressure changes due to freestream Mach number 
increase cannot be compensated by aft shock wave movement, thus resulting in a less rapid (smaller derivative value 
of dCLfdM _) rise in lift. By Mach 0.96 the entire upper surface shock wave is essentially at the trailing edge and 

the lift decreases, as indicated by the negative value of 3CLjdM m . As can be seen on the figure, the effects of this 

shock wave movement are captured by the variations in the sensitivity derivatives. Also, notice on Fig. 5 for 
supersonic freestream Mach numbers that the sensitivities are considerably lower. Additional results^* show that 
the derivatives of the total lift coefficient exhibit their largest change with respect to M a , T, C, a, XLj, XT t 
followed by T 2 , T 3 , L, Y T , T 4 , and T Jf indicating that a hierarchy of dominance exists among the design variables 
for the current wing configuration. Finally, again there is good agreement between the results obtained by the quasi- 
analytical method and the finite difference approach. 

Conclusion 

Based upon the above results, it is concluded that the quasi-analytical method is a feasible and efficient 
approach for accurately obtaining transonic aerodynamic sensitivity coefficients in three dimensions. In addition, 
use of the symbolic manipulation packages to carry out the symbolic evaluation of the elements of the sensitivity 
equations is crucial in this type of sensitivity study. The results obtained from the quasi-analytical method are almost 
identical to those obtained by the finite difference approach. Furthermore, the study indicates that: . 

(1) obtaining the sensitivity derivatives using the quasi-analytical approach and an iterative conjugate 
gradient method appears to be more efficient than computing the derivatives by the finite difference 
method, especially when the number of design variables is large, and 

(2) the quasi-analytical method shows promise with regard to analysis-sensitivity methodologies applied 
to large aerodynamic systems. 
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(a) Chordwise pressure coefficient 



(c) Sensitivity to angle of attack 



(e) Sensitivity to camber 



(b) Sensitivity to M m . 



(d) Sensitivity to thickness 
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9- 1 — Pressure coefficient and pressure coefficient sensitivity derivatives at 56% semispan, 
0.84, a = 3°. Values in lower left-hand corners are Cl for (a). For(b)-(f), they are the section 

computed by FD and QA methods. 
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Fig. 2 - Lift coefficient sensitivity derivatives, M m = 0.84, a = 3°; Values in lower left-hand 
corners are the wing computed by FD and QA methods. 
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Fig. 3 ~ Sensitivity to nonbasic design variables. Values in lower left-hand corners are the wing 
^ - computed by FD and QA methods. 








(a) Thickness increased to 0.065 chord 



(b) XLj moved aft 0.1 root chord 

Fig. 4 — Actual versus predicted pressure coefficients at 56% semispan. Values in lower 
left-hand corners are the actual and predicted section Cl. 
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Fig. 5 - Sensitivity of CL for 0.8 < M <12 
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NOMENCLATURE 

Chord length 
Maximum camber 
Section lift coefficient 
Wing lift coefficient 
Pressure coefficient 

Young's modulus of elasticity (non-dimensionalized by 10 7 Psi) 
Grid stretching factors 

Leading and trailing edge index in the x direction 

Location of maximum camber 

Freestream Mach number 

Residual of the aerodynamic equation 

Twist angle at the tip 

Maximum thickness 

Twist angle at a given section 

Thickness of the plate(non-dimensionalized by the chord) 
Residual of the structural equation 
Residual for the U equation 
Freestream velocity(also Uj n f) 

Perturbation velocity Cartesian components 
Vector of design variables 
Cartesian coordinates directions 
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V 

PoO 


Angle of attack 

Structural deflections non-dimensionalized by the chord 

Small perturbation velocity potential function 

Ratio of specific heats 

Circulation at a given station along the wing 

Poisson's ratio 

Freest ream density 


Subscripts 

00 

A 

b 

i, j>k 

ii, jj,kk 

iii, jjj,kkk 
LE 

I 

P 

root 
- S 
TE 
tip 
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Free stream condition 

Obtained from aerodynamic variables alone 

Body 

Grid point 

Counters for the residual dependencies 

Counters for the selected deflections and loads for the coupling 

Leading edge 

Lower side of the wing 

Pressure 

At the root 

Obtained from the structural variables alone 
Trailing edge 
At the wing tip 
Upper side of the wing 
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INTRODUCTION 

In the transonic regime, due to the non-linearity of the governing flow 
equations, the determination of optimum aerodynamic loads is one of the main 
difficulties facing the aircraft designer. Since most present day commercial aircraft 
operate transonically, computational methods which use optimization techniques are 
being developed to improve current designs. However, in order for these advanced 
computational codes to become more useful as design tools, it is necessary to develop 
methods for the computation of the sensitivity of the different parameters, such as 
aerodynamic forces or structural deflections, to the different design variables. With a 
sensitivity derivative being defined as a system response of interest with respect to a 
given independent design variable, it is desirable that such sensitivity coefficients be 
easily obtained. 

In the past, sensitivity methodology has been used in structural design 1 and 
optimization programs 2 and in some aerodynamic studies . 3 8 However, the 
predominant contributor to cost and computational time in the optimization procedure 
has always been the calculation of sensitivity derivatives. Hence, efficient numerical 
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methods for computing such derivatives are needed for the integration of advanced 
computational codes into systematic design methodology, where the computational 
cost of a single flow analysis can be extremely high, particularly in three dimensions. 

Consequently, the primary objective of this research is to investigate the 
concept that it is possible to use similar, perhaps identical, incremental iterative, 
solution approaches to efficiently couple for three dimensional transonic flow an 
aerodynamic solution for the pressure distribution with a structural solution for the 
corresponding deflections and to simultaneously use the same solution algorithms and 
the quasi-analytical method to obtain the aerodynamic as well as the structural 
discipline sensitivity derivatives for the fully coupled system with the input coefficients 
necessary to determine system sensitivities. Since the entire method is complex and 
requires an efficient flowfield as well as structural solver and since the present study is 
essentially proof-of-concept, it was decided for the present work to base the 
aerodynamics on the transonic small perturbation potential equation and the structural 
solver on the small deflection plate equation. Because of their simplicity, these 
equations are practical tools for the present proof-of-concept study where rapid 
solutions are essential. Previous experience with this approach has indicated that it is 
robust and reasonably accurate for engineering purposes. Finally, in order for an 
optimization process to be accurate, it must take into account the system sensitivity 
derivatives in which the effects of each discipline on the other is considered. Thus, the 
solver also computes the coupling derivatives relevant to the calculation of the system 
sensitivity derivatives. 

Currently, one conceptually simple method for computing sensitivity derivatives 
is the method of "brute force" finite differencing. Here, a design variable is perturbed 
from its previous value, a new complete solution is obtained, and the differences 
between the new and old solutions are used to obtain the sensitivity coefficients. This 
direct technique has the disadvantage of being very computer intensive, especially if the 
governing equations are expensive to solve and the number of design variables is large, 
and the resultant values are often very sensitive to the magnitude of the design variable 
perturbation. As a less costly alternative, sensitivity derivatives can, in principle, be 
computed by direct differentiation of the governing equations. In the case where the 
continuous governing equations are differentiated prior to their numerical 
discretization, the method is known as the "continuum" or the analytical approach 3 . On 
the other hand, if the governing equations are differentiated after their discretization, 
the method is known as the "discrete" or the "quasi-analytical" approach. 

Investigations concerning the feasibility of the quasi-analytical approach for the 
computation of the aerodynamic sensitivity derivatives have been undertaken by many 
researchers 4 - 5 ’ 6 and several methods have proven to be very successful. However, the 
differentiation of the governing discretized equations results in very large systems of 
algebraic linear sensitivity equations which must be solved to obtain the derivatives of 
interest. The application of a direct solver method to such a system requires extensive 
computer storage which for practical three dimensional problems is beyond the capacity 
of modem supercomputers. Moreover, the sensitivity matrix, sparse in nature, is 
generally very ill conditioned (or not diagonally dominant) and the convergence by the 
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use of standard iterative techniques is very slow. To avoid these problems, it is 
necessary to develop other iterative solution algorithms of the sensitivity equations. 
One possibility is the incremental iterative technique 4 - 7 which allows the iterative 
calculation of the sensitivity derivatives using algorithms similar to those applied to the 
flowfield.. 

The incremental iterative technique can be applied through a point semi-implicit 
algorithm to solve for the flowfield, structural deflections, and their respective 
sensitivities with respect to the different design variables simultaneously. However, 
these results are only discipline specific. To obtain a trully optimized solution the effect 
of one discipline on the other 8 needs to be considered. In other words, system as well 
as discipline derivatives need to be determined. Consequently, a second objective of 
the current work is to not only compute the coefficients needed for the system 
sensitivity equations but to also investigate the number of system sensitivities needed 
and methods for computing them. 

THEORY 

Flowfield Model 

The equations governing transonic flow are highly nonlinear and range from the 
Navier-Stokes equations to the small perturbation potential equation. Since this 
research is a proof-of-concept investigation, the flow modeling is the simplest possible, 
e.g. the non-conservative transonic small perturbation equation: 


(1- Mi - (y +1)M]> 0 0-a) 

where 

4>x=u (lb) 

4>y=v (l.c) 

<t> z =w . (Id) 


As shown in Fig. 1, the selected geometry is a rectangular wing with the z axis in the 
spanwise direction. 



Fig. 1 Geometry Setup 


At the wing, the boundary conditions are the inviscid surface boundary 
conditions for tangent flow: 
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/ 1 1 ^ . 5 y u . dy u . 

v - (U w + u) + w 


( 2 ) 


dx dz 

where y u , is a function of x and z and certain design variables, such as the angle of 
attack. In the wake, the Kutta condition along the wing semispan yields: 

r= A<}>, x TE < x <00 (3) 

while at the farfield, the boundary condition is: 

<t>oo=0 (4) 

At the downstream boundary, the Trefftz boundary condition can be approximated by: 

At x=oo : 4> x =0 (5) 

Further, the wing symmetry condition is expressed by: 

w(z=0)= 0 (6) 

The finite differencing of Eq. (1), requires the use of a residual R written in 
functional notation at the point i,j,k as: 

R.jjc-Riiuttii/w. XD) (7) 

Since the structural deflections are included as the boundary conditions, and are not 

treated as dependent design variables in the above equation, Eq. (7) should be 
considered a discipline equation. 

After taking the total differential of Eq. (7) with respect to a design variable 
XD, the sensitivity equation is obtained: 


4ft i.j.k 
dXD 


ii.jj.kk 


ii-ii-tt ] + 
dXD J 


i. j.lc 1 

<7XD 


= 0 


( 8 ) 


In this equation lies the essence of the quasi-analytical formulation in which the 
discretized governing equations are differentiated. Here, ^ is <{>(x(ii), y(jj)> z(kk), 
XD, 5); and the system matrix dR/d<j) is sparse, or non zero at certain points only 
(mostly the ones neighboring i,j,k). In this equation, the vector of deflections {5}, even 
though not explicitly shown, is considered to be a vector of independent variables. 
Near the boundaries, Eq. (8) has been reformulated to include the flowfield boundary 
conditions. The flowfield sensitivity derivatives d§/dXD that are obtained from solving 
Eq. (8) above can be used to calculate pressure sensitivities dCp/dXD which in turn can 
be used to calculate the sensitivities of the section and wing lift to the design variables. 


Structural Modeling 

The structural problem is modeled by representing the wing by an equivalent 
flat plate with dimensions almost coincident with those of the wing. The equation 
describing the plate deflections is: 9 

DV 4 5-Ap=0 (9) 

which assumes a thin plate and small deflections. Here, Ap is the loading due to the 
difference in pressures between the upper and lower surface: 

Ap = J-/O.ULAC, and AC„ = C,,-C W (10) 
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and D, the flexural rigidity (or equivalent bending stiffness), is given by: 


- ill] 

12(1- v 2 ) 

This model, while simple, will yield both bending and twisting effects. 

The boundary conditions 10 for Eq. (9) involve both fixed and the free edges. 
The root is the only fixed edge and there the boundary condition is: 

At z=Q 8= 0 (12. a) 


At the tip, the boundary conditions are written: 


(12. b) 


■tip • n 3 + (^ 


dz dx 


(13 a) 


v ^r r = 0 (13-b) 

Eq. (13. a) combines the no twisting moment and no shearing force conditions at a free 
edge while Eq. (13 .b) states that the bending moment along the edge is zero. The other 
two free edges are the leading and trailing edges, and the boundary conditions for those 
are written: 


x = X LE > X TC 


^+(2-v)-^V=0 (14. a) 

dxdz 


f4- + 0 (,4.b) 

dx 2 dz 2 

Hence, this system of equations establishes a well defined boundary value problem that 
can be solved by finite differencing. 

The residual for Eq. (9) can be expressed as: 

T lk =T ik (8 likk ,XD) (15) 

Again, this equation is discipline specific since 8jj ^^(x^i), z(kk), <j> jiikkk , XD); and 
^iii.kkk 1S vector of potentials on the upper and lower side of the wing that are 
related to the calculation of loads. This vector is considered to be composed of 
independent variables. Unlike the flowfield case, which is three dimensional, the 
deflection field is a two dimensional variable. 

After taking the total derivative of Eq. (15) with respect to a design variable the 
structural sensitivity equation is obtained: 


dXD 


dXD 


= 0 ( 16 ) 
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In this case, also, the system matrix is sparse. Like the flowfield case, Eq. (16) must 
take into account the appropriate conditions at the boundaries. 


Cou pling 

As a result of aerodynamic loads, the equivalent plate representing the wing will 
deflect; and such deflections will perturb through bending and twisting of the wing the 
section angles of attack and camber line shapes. These deflections in turn will induce 
different load distributions, and the two processes must be interacted until a converged 
solution is obtained. This interaction is the process of aerodynamic and structural 
coupling. 

The coupling between the structural and the flowfield solutions is achieved 
through the wing boundary conditions and is included by simply adding the structural 
deflections to the ordinates of the wing. Hence, after taking the derivatives with 
respect to the x and z coordinates, the boundary conditions equations are modified by. 


- ( il "•' ) + (17.a) 

dx V J A 

^ = f *?■•'. ) + £L (17.b) 

3t. I 3z J A 3z 

Note that this coupling is only carried out at the field variables level. In other words, 
for a linear case (much below the critical Mach number), a case in which the sensitivity 
matrix dR/dfy would not be influenced by the values of (J), this aeroelastic coupling 
would only slightly affect the aerodynamic and the structural sensitivity derivatives. 
Thus, the coupling is said to be achieved for the sensitivity derivatives at the zero order 
only. 


System Sensitivity 

As mentioned, for an optimization process to be accurate, it must take into 
account the system sensitivity derivatives in which the coupling between the disciplines 
is included. Thus, the calculation of interdisciplinary sensitivities such as the sensitivity 
of the pressure distribution to the thickness of the plate or that of the tip deflections 
with respect to the camber at the tip are needed. In general, the set of equations 
governing the entire coupled system can be written as: 9 

A((XD, 5), <J>)= 0 (18) 

S((XD, (j>), S)= 0 (19) 

where Eq. (18) represents the aerodynamics and Eq. (19) is for the structures. For the 
system analysis <t> can be replaced by AC p since it is the variable involved in the 
aerodynamic coupling. The vectors grouped in the inner parentheses are the input, 
while the vectors of unknowns (output) are listed last. The purpose of the analysis is to 
find the total derivatives dY/dXD of the output vector with respect to the different 
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design variables. According to the implicit functions theorem, the equations above can 
be written as 10 : 

AC P = AC Pa (XD, 5) (20) 

5= 5 S (XD, ACp) (21) 

After considering Y=((ACp),{5}), taking its total derivatives with respect to XD, and 
rearranging the terms, the following system equation is obtained: 


L-J sa I JldXD j [ dXD J 

where J^yg is a Jacobian of the partial "coupling" sensitivity derivatives d&Cp/dd and 
Jsa is the Jacobian of d5/3ACp for selected points on the wing. For example, the i-th 
column of J^g comprises the partial derivatives with respect to the i-th displacement. 
The partial derivatives in the coupling matrix as well as the right hand side are, by 
definition, calculated using strictly discipline derivatives. Again, the quasi-analytical 
approach is used. Equations (7) and (15) are rewritten as: 

5),4. iijiw ) (23) 

T Ul =T Ul ((XD,ACp).5 ii , ldt ) (24) 

where 8 is considered an independent variable for Eq. (23) and ACp is considered 
independent for Eq. (24). This approach is valid since it is discipline specific. 
Differentiating Eq. (23) with respect to a given deflection and Eq. (24) with respect to 
a given ACp on the wing yields the system of linear coupling sensitivities equations. 


X i.j.fc _ 

d 5 


<¥ ii.jj.v 


ii.jj.kk 

d5 


i.j.fc 1 _ 


d A C 


<?A C 


<2 A C 


which when solved yields the coupling sensitivity derivatives, i.e. the elements of the 
J^ and J SA matrices, necessary for the calculation of the system sensitivities via 
Eq. (22). 


DESIGN VARIABLES 


Design variables are classified into two groups, the aerodyanmic variables 
termed XDA and the structural variables called XDS. One variable (M^) is common 
to both vectors. A design variable is termed to be aerodynamic or structural depending 
in which expression of the discipline residuals it appears. For example, the angle of 
attack would be an aerodynamic variable while the plate thickness, which only appears 
in the deflection equation would be a structural variable. However, all the design 


variables used are basic variables in that they are uncoupled and independent. For the 
current problem the vector of design variables consists of twelve variables and ts given 

by: 

XD= (XD1,..., XD12) (27) 

These design variables can be classifi ed into three groups. 

(a) Freestream design variables: These include the freestream Mach number and the 
angle of attack. The Mach number enters the formulation through Eq. (l) while a 
appears in the boundary conditions in Eq. (2). 

(b) Cross section design variables: These include the variables that define the airfoil. 
For the present study only NACA four-digit airfoils are considered. Thus the relevant 
design variables are maximum thickness, ma xi mum camber, and location of the 
maximum camber, at both the root and tip. Another variable defining each spanwise 
section would be the geometric twist, usually defi ned. m te rms of the relative twist of 
the wing tip to that of the root. The airfoil sections as well as the aerodynamic twist at 
a given span station are obtained by linear lofting between the root and tip the values 
for TH, C, and LC, each expressed as a fraction of the chord, i.e.: 

TH = TH root +z/z ti p(TH ti p-TH root ) (30) 

C = C root +z/z t ip (Mip- C roo t) 

LC — LC roo t+z/z t jp (LC t ip- LC root ) ( ) 

It should be noted that this formulation is not a point by point lofting in which the 
vertical coordinate is interpolated linearly from root to tip. Nevertheless, this approach 
was chosen to simplify the analytical derivations as well as the coding. The section 
twist is also obtained by linear interpolation between the wing root and the wing tip: 

tw= z/z t j p T t j p ( 33 ) 

When taking the derivative of y u , with respect to x, the following is obtained. 


dy 


X ) A 


dy c . dy 


dx 


TH 


dx 


-a - tw 


(34) 


_ dy * , dy 


TH 


Tu p x 


dz 


dz 


dz 


(35) 


A ~ ~ ~ “ li P 

With this formulation, the vector of the aerodynamic design variables can be written: 
XDA= (a, TH roo t, TH^jp, C roo t> C^p, LC r0 ot, ^Cfip> Tfip> ^°o) 

(c) The structural variables: These include the parameters involved in the plate 

deflection equation. The first, M*,, comes from the dynamic pressure term, second is 

the thickness of the plate t, followed by Poisson's ratio v, and Young's modulus of 
elasticity E. In the present study, the dynamic pressure is calculated using the sea level 
conditions. Thus, the vector of structural design variables is. 

XDS=(M„, t, v, E) (37) 

These two vectors are combined to form a single vector of design variables: 

XD= (a, TH r00 (, TH t jp, C r00 (, C^jp, LC roo t > LC^jp, T t jp, M^ t, v, E) (38) 
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DISCRETIZATION APPROACH AND NUMERICAL PROCEDURE 


Aerodynamic Analysis 

As previously stated, the aerodynamic analysis is based on the transonic small 
perturbation potential formulation in Eq. (1), formulated using a Cartesian grid and a 
finite computational domain. Hence, the transformation utilized maps the infinite 
physical domain into a finite computational grid. 

In the present formulation, the infinite physical plane is transformed via tangent 
functions into the finite computational space shown in Fig. 2. Thus, the i=l, i=IM, j=T, 
j=JM, and k=KM planes physically correspond to infinity and k=I is the wing symetry 
plane. Further, the wing is located between two grid lines, JB and JB-1 . 



Fig. 2 Computational and Physical Domains 


Since there is a potential jump across the trailing edge cut, which extends to 
downstream infinity, and since the jump depends upon the trailing edge potentials, the 
sensitivity matrix, dR/d<f>, while banded, contains elements far from the central band. 
Consequently, the rapid and efficient solution of the sensitivity equation by direct 
methods is difficult. However, the sensitivity equation can be solved by the same 
iterative method as the flowfield, by the introduction of a new residual S;j k 
corresponding to Eq. (8). 

For the sensitivity portion of this analysis, the residual expressions, Rg k , are 
differentiated analytically with respect to the flowfield variable <j>. Similarly, the right 
hand side of the sensitivity equations is determined by analytical differentiation of the 
residuals with respect to each design variable. Unfortunately, the size of the sensitivity 
matrix is tremendous for fine grids and storing such a matrix is beyond the capability of 
many computers. In the present study, the storage problem is solved by the use of a 
more efficient solver, namely the incremental iterative technique. The additional 





available memory space allows the use of finer grids and the inclusion of other 
disciplines as well. 

Once the potentials are obtained, the C p 's are obtained through the small 
perturbation relation: 

C p = -2<j> x (62) 

which requires the extrapolation (first order) of the <J> values above and below the plane 
of the wing. The section lift coefficient is then computed directly from the circulation 
around every airfoil section: 

C| = 2T k = 2(<j> ITE ro k — ^ITE.JB-I.k ) (63) 

which gives faster and more accurate results than those obtained by integration of the 
C p 's difference between the lower and upper surfaces. The wing lift coefficient C L is 
calculated by integration over the span. The corresponding sensitivity derivatives are 
then determined from: 


ac 


■1 = 


ar, 


(64) 


V axD ; k axD 
and aC L /aXD can be calculated by numerical integration of the above coefficient along 
the span. 


Structural Analysis 

Since Eq. (9) is a fourth order partial differential equation, its solution can be 
significantly simplified by splitting it into two equations to be solved simultaneously: 11 


In non-dimensional form, Eq. (65) is: 


where 


and 


ko = 


6yp 


1 > 
|T3 

II 

o 

(65) 

D 

u = 0 

(66) 

, = 0 

(67) 

- v’)“- 

(68) 

0 

(69) 


After splitting the governing equations, the next step should be to split the eight 
boundary conditions written in Eqs. (12)-(14) into four for the u's and four for the 
deflections. However, this splitting must be carried out so that the solution scheme 
does not become unstable. When applying the Laplacian operator to Eq. (12. b) one 
obtains: 
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At z=0: 


d u 

d z 


= 0 


(70) 


For the free edges, the situation is more complicated. The first step would be to take, 
the partial derivative of u in the z direction. When combining the result of this 
differentiation with Eq. (13. a) the appropriate boundary condition becomes: 


At z=z ,„ 


du 


= (V- 1) 


d 3 5 


(71) 


1lp dz ' dzdx 2 

When finite differenced, this formulation is accurate and has better stability 
characteristics than the formulation with the third order partial derivative with respect 
to z. Similarly, at the leading and trailing edges of the plate the boundary condition for 
u can be expressed: 


X X LE> X TE : 


du 


= (V - 1) 


d 3 8 


(72) 


dx ' dx d z 2 

Since the flowfield solver uses a finite differencing technique with a given grid 
the same technique and the same grid are used to obtain the structural deflections, 
which simplifies the aerodynamic structural coupling. Hence, both the flowfield and 
structural solutions can be calculated within the same loops, which is computationally 
efficient. Consequently, the structural part uses the same grid metrics. Further, the 
field variable u is first obtained and then used as an input to solve for the deflection at 
the same point, thus enhancing convergence and stability. 

The boundary conditions stated in Eqs. (12)-(14) should be applied at the exact 
boundaries of the wing which do not coincide with an exact grid point in the 
computational domain because the leading, trailing, and wing tip edges are located 
between grid lines. Hence, a Taylor series development should be used at all the 
boundaries except the root, where the boundary coincides with a gridline. This 
development would involve higher order partial derivatives which when finite 
differenced would yield extremely complicated expressions. To avoid that problem, the 
size of the equivalent plate and the grid were chosen such that the free edge boundaries 
of the equivalent plate are very close but not coincident to the wing leading, trailing, 
and tip edges. The boundary conditions corresponding to Eq. (71) and (72) can cause 
numerical divergence and a possible solution is to simplify them so that numerical 
stability can be created. Fortunately, the variable u, physically corresponds to the 
second derivative in one direction at a given edge. For example, at the wing tip, if it is 
assumed that the loading along the wing is only a distributed loading without 
concentrated loads or moments to cause discontinuities in the curvature of the plate, 
the assumption that the second partial derivative with respect to z is constant is 
acceptable. Similarly, at the tip, the "curvature in the x-direction" or the second partial 
with respect to x will also be constant in the absence of concentrated loads and 
moments. Hence, the approximation that u is constant along a free tip is reasonable. 
In addition, this condition does not have a destabilizing influence on the algorithm. 
Hence, the partials with respect to z and x respectively for Eq. (71) and (72) are 
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assumed to be equal to zero. Also, the boundary conditions for the structural 
deflections equation (Equ. 66), stated in Equ. (12 a), (13.b), and (H.b), are finite 
differenced and incorporated in the corresponding residual expressions. 

The second part of the structural analysis is the structural sensitivity analysis 
with respect to the four components of XDS, the vector of structural design variables. 
The approach used is the same as the one used for the deflections. In other words, the 
sensitivity equation is also divided into two components. Hence, when applying the 
quasi-analytical approach to Eq. (75) and (78) the following equations are obtained: 


STU 


i,k 


STU 


i ,k 


d u 


it , kk 


g u ■■ .a 1 + 
d XDS | 


8TU JJ_1 = 0 
d XDS 


( 88 ) 


ST 


i.k 


u 56 


^ u .kk 1 + 

a xds J 


a xds 


(89) 


Here, it should be noted that 


dT,, 


axDs 


is 


axDs 


("•!) 


which shows that, as in the 


deflection field solution, the output variable of the system of Eq. (88) is used as an 
input to Eq. (89). At the boundaries, Eqs. (88) and (89) must take into account the 
appropriate structural boundary conditions. 


Aeroelastic Coupling and System Sensitivity Analyses 

Aerodynamic-structural coupling can be carried out at two levels; defined here 
to be zero and first order. The zero order coupling corresponds to an updating of the 
aerodynamic boundary conditions each time, after the structural deflections are 
calculated and vice versa. However, sensitivities are computed as discipline 
sensitivities and do not directly include the complete effects of aerodynamic-structural 
coupling. On the other hand, the first order coupling is defined to mean that the effect 
of the structure on the flowfield and vice versa is taken into account not only at the 
flowfield-deflection level but also at the sensitivity level. For example, for the zero 
order coupling the structural deflections affect the aerodynamic sensitivity denvat.ves 
through the spanwise flow component 4> z in SRJdXDA while the first order coupling 
also affects that expression through a coupling term dfyJdS. This term is called a 
coupling sensitivity. In this second case, the deflections are not considered constant in 
the aerodynamic residual expressions (Eq. 18), as in the discipline specific analysis, and 
are considered as design variables. Likewise, in first order coupling, the potentials 
related to the C calculation along the wing are treated as design variables for Eq. (19). 

The terms that affect coupling the most are those that appear directly in the 
residuals expressions. These are the deflections, since they enter directly in the 
expression of the boundary conditions for the aerodynamic residuals, and the loads 
AC p , which appear in the expressions of the structural residuals. However, as shown in 
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Eq. (22), the coupling derivatives, 5ACp/c)5 and <35/dACp, are the essential components 
of the system sensitivity matrix and are used to obtain the system sensitivity derivatives. 
The equations for determining those coupling derivatives are presented in Eq. (25) and 
(26). However, frequently not all the deflections or loads can be used in the system 
matrix since such inclusions would often require extensive memory storage and CPU 
times that are unrealistic. Hence, the choice of which loads and deflections to include 
in the system sensitivity equation is subject to judgment and experimentation. 
However, the more coupling variables are included, the more accurate the system 
sensitivity derivatives should be. 


Numerical Approach 

The sensitivity matrix, associated with the linear sensitivity equations, as well as 
the matrix resulting from the finite differencing of the flowfield and structural solutions, 
are generally very sparse and ill conditioned, or not diagonally dominant. Thus, the 
solution of the corresponding linear equations by standard direct solvers is memory 
inefficient and iterative methods should be considered. 6 - 7 ’* In addition, since the non- 
linear flowfield equations must be solved iteratively, the use of a similar iterative 
scheme to obtain the sensitivities would seem to be appropriate. 

A possible scheme is the incremental iterative technique, 4 - 8 which has shown to 
have better convergence characteristics in many cases than the standard iterative 
techniques. This method comes from a formulation in which a system of algebraic 
equations has the general form: 8 

[A]{z*} + { B ) = {o} ( 90 ) 

where {Z*j, the solution vector, is obtained by the two step formulation: 

-[a]{a-z}= [A]{Z"} + {B] (91) 

{ z "*' } = (z"} + {a"z} (92) 

Here, n is the iteration index and [z] is a convenient approximation of [Z], generally 

chosen to enhance the diagonal dominance and, thus, the convergence characteristics of 
the system. 

The above formulation, when applied to sensitivity equations, still requires the 
storage of a relatively large sensitivity matrix. However, the use of a point algorithm to 
obtain the increments avoids that problem since it only requires the elements of the 
matrix relevant to the calculation of the increment at point i,j,k. Obviously, such an 
approach has the possible disadvantage of slower convergence. Nevertheless, since the 
sensitivity equations are linear, their convergence should be faster than that of the 
nonlinear flowfield. Unfortunately, the structural equations tend to behave like the 
nonlinear flowfield equations in terms of convergence. 

An example of such a point algorithm is the semi-implicit ZEBRA scheme 12 
which mimics point successive over relaxation (SOR). The algorithm marches in the 


streamwise (I) direction solving by spanwise planes. In each plane, the points where 
j+k are odd are denoted black and the ones where j+k are even are denoted white (Fig. 
3). Each plane is solved by a two-pass sweep in which new black values are obtained 
First, followed by the white ones. Convergence is thus accelerated because calculations 
at the white points use updated values at the black points. Because of its uncoupled 
formulation, this method is suitable for sequential, vector, and parallel machines. 



In the ZEBRA algorithm, because of the point semi-implicitness, the matrix [^] 

is reduced to a scalar B. Hence, the incremental changes in the unknowns can be found 
in the following form for the aerodynamic potential, for example: 

A4>.. fc =^- + DMP (93) 

,J - B 

where DMP is a damping term added for transonic stability. The same type of formulas 
can also be used to calculate the increments for the aerodynamic sensitivity field 
variables, structural deflections, structural sensitivity derivatives and coupling 
derivatives field variables^. The algorithm used is schematically described in the 
following figure. 
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Iteration Loop 
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•Use results to update aerodynamic and structural boundary conditions 
►Iterate until convergence 


•Solve System Sensitivity Equations for 


dzNCp d5 

and 

dXD dXD 


Fig. 4 Integrated Solution Approach 


RESULTS AND DISCUSSION 


Case Studied 

The wing configuration considered in this study has a rectangular planform of 
constant unit chord. The geometric and structural design parameters describing the 
wing as well as the nominal freestream conditions are listed in Table 1. It should be 
noted that the root and tip airfoil sections are NACA 2406 and NACA 1706 
respectively. 
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Table! Wing Data 


Aspect Ratio 

3.17 

a 

2.0° 

^ip 

1.58 

c 

1.0 

T ..p 

-1.0° 

TH rooI 

0.06 

TH lip 

0.06 

c 

v “'rool 

0.02 

Cup 

0.01 

LCfOOt 

0.4 

LC tip 

0.7 

M„ 

0.82 

E 

1.0 

v 

0.33 



The case presented here used a grid of 97 x 16 x 16 for the flowfield and 49 x 
10 for the structural deflections. The freestream Mach number is supercritical and a 
shock wave is present on the inboard sections of the upper surface of the wing. 
However, the shock wave disappears on the outb oard sections due to three dimensional 
effects. Thus, this Mach number is interesting because it locally includes both the 
subcritical and supercritical behavior of the flowfield and the corresponding 
sensitivities. In all cases, to speed up convergence, a coarse-medium-fine grid 
sequence halving in the x-direction was used in computing the analysis information. 
Results were computed for equivalent plate thicknesses of five and two percent but 
only the two percent results were shown in this paper. It should be noted that a one 
percent thick case, while attempted, turned out to be aeroelastically divergent. For the 
coupling variables needed to determine the system sensitivities, five of the ten spanwise 
stations were selected each involving twenty five of the forty-nine possible points. It is 
believed, since it is numerically difficult to include every point used in the fine grid, that 
the deflections and loads selected for the coupling system coefficients in the sensitivites 
will be representative. However, this choice is under investigation and will be further 
discussed in the final paper. 

Fig. 5 shows the pressure distribution at six of the ten spanwise stations. The 
upper surface shows a shock wave at approximately the x=0.5 location in the sections 
near the root. The airfoil section, being non-constant from root to tip, is also drawn on 
the same diagram and the angle of attack as well as the geometric twist are taken into 
account when plotting the geometry. The final deflected shape of the airfoil due to 
aeroelastic coupling is drawn in dashed lines but not to scale. The critical pressure 
coefficient level C p \ is also shown and comparison of C p * with the pressures shows 
that the shock wave weakens progressively when approaching the tip, which is 
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obviously subcritical. At the tip, the pressure distribution is typical of a subcritical afl 
cambered NACA section. One should note that due to the change in airfoil sections 
from root to tip the wing has some inherent aerodynamic twist. However, unlike the 
thick plate case 16 , the lower surface Cp curve at the tip section, goes above the upper 
one causing the aerodynamic load at the leading edge to be negative. 

The results for the discipline sensitivity derivatives are shown in Fig. 6-14 and 
Fig. 16-19 while the ones for the system sensitivity derivatives are shown in Fig. 20 
through 3 1 . Pertinent portions will be selected and discussed in detail in the final 
paper. 

Fig. 15 shows the structural deflections at different span stations. Notice that 
if a line is drawn from the leading to the trailing edge of the plate at each section, this 
line would form an "angle of attack" with the x-axis which would be an induced twist 
due to structural deflections. Further, even though the amplitudes are extremely small, 
bending exists in the sections toward the wing tip. This "cambering" effect due to 
chordwise bending is more pronounced as the tip is approached. In fact, the chordwise 
section of the equivalent plate near the tip looks as a camber line that could cause an 
increase in lift and could become a dominating component of the tip aerodynamics. 
Further, the maximum of the structurally induced camber is a little bit aft of center. 
Note that the spanwise edge of the equivalent plate is loaded due to the ACp there, 
even though no concentrated loads or bending moments exist at the edges of the 
equivalent plate. If the spanwise edge of the equivalent plate actually corresponded to 
the wing tip, it would not be loaded and the cambering effect would be attenuated at 
the tip. Again, pertinent features of the structural sensitivity derivatives shown on 
Figs. 16-19 will be discussed in the final paper. 

For clarity and length reasons, the system sensitivity derivatives plots are only 
shown for three stations, in Fig. 20 to 3 1 . For the sensitivity of the loads with respect 
to the design variables the system and discipline curves almost agreed. However, a 
discrepancy was noticed at the leading and trailing edge locations for the sensitivity of 
the loads. Moreover, when compared to a case where thirteen stations chordwise were 
chosen for eight spanwise stations differences were found for the structural system 
sensitivities with respect to the aerodynamic design variables and for the loads system 
sensitivities with respect to the structural design variables. This difference is currently 
under study and will be discussed in the final paper. However, in all cases the system 
structural sensitivities, d5/dXD, often differed in magnitude and, more importantly in 
sign from the discipline sensitivity derivatives, 05/3XD. The origins and significance of 
this behavior will be discussed in the final paper. 
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Fig. 6 Pressure Coefficient 
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Fig. 7 Pressure Coefficient Sensitivity to the Root Thickness. 
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Fig. 8 Pressure Coefficient 
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Fig. 9. Pressure Coefficient 
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Fig. 10 Pressure Coefficient Sensitivity to the Tip Maximum Camber. 
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Fig. 1 1 Pressure Coefficient 
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Fig. .13 Pressure Coefficient Sensitivity to the Twist at the Tip. 
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Fig. 14 Pressure Coefficient Sensitivity to the Mach Number. 
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Fig. 22 System Sensitivity Results for the Tip Thickness. 
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Fig. 23 System Sensitivity Results for the Root Camber. 
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Fig. 24 System Sensitivity Results 
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Fig. 27 System Sensitivity Results for the Tip Twist. 
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Fig. 30 System Sensitivity Results 
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Summary 

The system sensitivity analysis shows that the deflections at the tip as well as 
the loads are going to decrease with an increase in a. Likewise, an increase in the 
values of TH rool , TH lip , C rool , C tip , LC rool , LC tip , T lip , and will cause an unloading of 
the tip associated with a decrease in the deflections. Further, the structural variables t, 
v, and E will cause an increase in the tip loading as well as the associated deflections. 
The structurally induced camber is essential to the interpretation of these a priori 
unexpected results. 


CONCLUSION 

Based on the results presented, the use of the incremental-iterative technique 
through the semi-implicit ZEBRA scheme to calculate the sensitivity derivatives 
obtained from the quasi-analytical formulation has proven to be successful and very 
computationally efficient. A large memory space for the storage of the sensitivity 
matrix is not needed anymore and the sensitivity derivatives can be calculated at the 
same time as the flowfield instead of using a converged flowfield solution as an input to 
a sparse matrix solver . 6 

The saved computational resources can thus be used for finer grids, more 
design variables, and additional disciplines. Hence, a coupling of the aerodynamic 
solver with a structural one and its sensitivities has been undertaken. This static 
aeroelastic coupling is very efficient since the structural calculation and resultant 
structural sensitivities and coupling sensitivities are computed at the same time as the 
flowfield. In addition, the use of finite differencing to solve for the structural 
deflections improves the efficiency of the scheme since no grid transformation is 
necessary. 

Because the system is multidisciplinary, the calculation of the system sensitivity 
derivatives takes into account the influence of one discipline on the other. This 
calculation relies on the "coupling" sensitivity derivatives that are not very easy to 
obtain computationally since their respective convergence (especially for the deflection 
with respect to load sensitivity) is slow. Results for the system sensitivities will be 
discussed in the final paper. 
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